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 toughest part. Now, we now have the better task of organizing what we did in a way that makes it general, hiding the small print while keeping the essential elements visible. In a way, we would like the optimizer to seem like a “magic box” that even users not aware of math modeling are capable of use to unravel their TSP problems intuitively.

4.1. TravelingSalesmanOptimizer design

Our optimizer class can have “core” methods, doing the majority of the work, and “superficial” methods, serving because the high-level interface of the category, which invoke the core methods underneath.

These are the steps that may lie on the core of the optimizer’s logic:

1. Create a Pyomo model out of a distance matrix. This is completed by the _create_model method, which principally 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 one essential difference between what we did and what we’re doing is that, now, the initial site is just not hard-coded as simply "hotel", but is assumed to be the positioning of the primary row in df_distances. In the final case, thus, the initial site is taken to be the primary one within the coordinates dataframe⁴ df_sites. This generalization allows the optimizer to unravel any instance.
2. (Try to) Solve the model. That is performed within the _optimize method inherited from BaseOptimizer, which returns True provided that an answer is found.
3. Extract the answer from the model and parse it in a way that is simple to interpret and use. This happens inside _store_solution_from_model, which is a technique that inspects the solved model and extracts the values of the choice variables, and the worth of the target function, to create the attributes tour_ and tour_distance_, respectively. This method gets invoked provided that an answer exists, so if no solution is found, the “solution attributes” tour_ and tour_distance_ never get created. The advantage of that is that the presence of those two “solution attributes”, after fitting, will inform the user of the existence of an answer. As a plus, the optimal values of each the variables and objective will be conveniently retrieved at any point, not necessarily in the meanwhile of fitting.

The last 2 steps — finding and extracting the answer — are wrapped contained in the last “core” method, _fit_to_distances.

“But wait” — you would possibly think — “Because the name implies, _fit_to_distances requires distances as input; is not our goal to unravel TSP problems using only coordinates, not distances?”. Yes, that is where the fit method matches in. We pass coordinates to it, and we make the most of GeoAnalyzer to construct the space matrix, which is then processed normally by _fit_to_distances. In this manner, if the user doesn’t need to collect the distances himself, he can delegate the duty by utilizing fit. If, nevertheless, he prefers to make use of custom data, he can assemble it in a df_distances and pass it to _fit_to_distances as a substitute.

4.2. TravelingSalesmanOptimizer implementation

Let’s follow the design outlined above to incrementally construct the optimizer. First, a minimalist version that just builds a model and solves it — with none solution parsing yet. Notice how the __repr__ method allows us to know the name and number of websites the optimizer comprises every 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 will be stated like: "Given a set of locations (and frequently 
their pair-wise distances), find the tour of minimal distance that 
traverses all of them exactly once and ends at the identical location 
it began from. For a derivation of the mathematical model, see [2].
Parameters
----------
name : str
Optional name to offer to a specific TSP instance
Attributes
----------
tour_ : list
List of locations sorted in visit order, obtained after fitting.
To avoid duplicity, the last site within the list is just not the initial 
one, however 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 needs to be picked because the "initial" one, doesn't matter which 
# really; by lack of higher criteria, take first site in dataframe 
# because 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 space 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, inside=pyo.Binary, 
doc="Whether to go from site i to site j (𝛿𝑖𝑗)")
model.rank_i = pyo.Var(model.sites_except_initial, 
inside=pyo.NonNegativeReals, 
bounds=(1, len(model.sites_except_initial)), 
doc=("Rank of every site to trace 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 have to 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 need to departure to precisely 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 every pair of non-initial sites (i, j), 
if site j is visited from site i, the rank of j have to be 
strictly greater than the rank of i.
"""
if i == j:  # if sites coincide, skip making 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:
""" Variety of locations to go to """
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 doesn’t contain any number of websites, because the representation string shows, or an internal model, and it’s in fact 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 space data, and if we don’t get a warning, it signifies that all of it went well. We will 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 space 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 within 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 within the optimal tour. If you happen to recall from their definition, they’re defined over all sites except the initial one⁵, and that’s why the hotel doesn’t appear in them. Easy, we could add the hotel with rank 0, and there we might have the reply to our problem. We don’t have to extract 𝛿ᵢⱼ, the choice variables for the individual arcs of the tour, to know by which order we must always visit the sites. Although that’s true, we’re still going to make use of the arc variables 𝛿ᵢⱼ to extract the precise sequence of stops from the solved model.

💡 Agile doesn’t have to be fragile

If our only aim were to unravel the TSP, without seeking to extend the model to encompass more details of our real-life problem, it might be enough to make use of the rank variables to extract the optimal tour. Nonetheless, because the TSP is just the initial prototype of what is going to change into a more sophisticated model, we’re higher off extracting the answer from the arc decision variables 𝛿ᵢⱼ, as they can 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 subsequent articles, selecting 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 go to some sites. Hence, if we extract the answer from 𝛿ᵢⱼ, our approach can be general and re-usable, irrespective of how complex the model we’re using.

The advantages of this approach will change into apparent in the next articles, where recent requirements are added, and thus additional variables are needed contained in the model. With the why covered, let’s jump into the how.

4.2.1 Extracting the optimal tour from the model

• We’ve got the variable 𝛿ᵢⱼ, indexed by possible arcs (i, j), where 𝛿ᵢⱼ=0 means the arc is just not chosen and 𝛿ᵢⱼ=1 means the arc is chosen.
• We wish a dataframe where the sites are within the index (as in our input df_sites), and where the stop number is indicated within the column "visit_order".
• We write a way 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 chosen arcs from 𝛿ᵢⱼ, which represents the tour. Done in _get_selected_arcs_from_model.
2. Convert the list of arcs (edges) into an inventory of stops (nodes). Done in _get_stops_order_list.
3. Convert the list of stops right into a dataframe with a consistent structure. Done in _get_tour_stops_dataframe.

As the chosen arcs are mixed (i.e., not in “traversing order”), getting an inventory of ordered stops is just not that straight-forward. To avoid convoluted code, we exploit the proven fact that the arcs represent a graph, and we use the graph object G_tour to traverse the tour nodes so as, arriving on 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 choice 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 chosen arcs from the choice 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 an inventory **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, because 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 take a look at again")

Let’s see what this recent method gives us:

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

The visit_order column indicates we must always go from the hotel to Notre Dame, then to the Louvre, and so forth, 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 now have the answer in a format easy to interpret and work with. However the sequence of stops is just not all we care about. The worth of the target function can be a vital metric to maintain track of, as it is the criterion guiding our decisions. For our particular case of the TSP, this implies getting the overall distance of the optimal tour.

4.2.2 Extracting the optimal objective from the model

In the identical manner that we didn’t use the rank variables to extract the sequence of stops because in additional complex models their values wouldn’t coincide with the tour stops, we won’t use the target function directly to acquire the overall distance of the tour, although, here too, each measures are equivalent. In additional complex models, the target function can even incorporate other targets, so this equivalence will not hold.

For now, we’ll keep it easy and create a non-public method, _get_tour_total_distance, which clearly indicates the intent. The main points of where this distance comes from are hidden, and can depend upon the actual targets that more advanced models care about. For now, the small print are easy: get the target 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 overall distance of the optimal tour"""
if self.is_fitted:
# as the target is an expression for the overall distance, 
distance_tour = self.model.objective()  # we just get its value
return distance_tour
print("Optimizer is just not fitted to any data, no optimal objective exists.")
return None

It could look superfluous now, however it’ll function a reminder to our future selves that there’s a design for grabbing objective values we’d higher 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 each the optimal tour and its distance are essential, let’s make the optimizer store them together every time the _fit_to_distances method gets called, and only when an optimal solution is found.

4.2.3 Storing the answer in attributes

Within the implementation of _fit_to_distances above, we just created a model and solved it, we didn’t do any parsing of the answer stored contained in the model. Now, we’ll modify _fit_to_distances in order that when the model solution succeeds, two recent attributes are created and made available with the 2 relevant parts of the answer: the tour_ and the tour_distance_. To make it easy, the tour_ attribute won’t return the dataframe we did earlier, it would return the list with ordered stops. The brand 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 space matrix 
provided in `df_distances`, after which optimizes it. 
If the model has an optimal solution, it's extracted, parsed and 
stored internally so it may possibly be retrieved.
Parameters
----------
df_distances : pd.DataFrame
Pandas dataframe where the indices and columns are the "cities" 
(or any site of interest) of the issue, and the cells of the 
dataframe contain the pair-wise distances between the cities, i.e.,
df_distances.at[i, j] comprises the space 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 an answer 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 space data and see how easy it’s to get the answer 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 will do even higher. To further increase the usability of this class, let’s allow the user to unravel the issue by only providing the dataframe of websites coordinates. As not everyone will have the option to gather a distance matrix for his or her sites of interest, the category can handle it and supply an approximate distance matrix. This was done above in section 3.2 with the GeoAnalyzer, here we just put it under the brand 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 salesperson" desires to visit, having the 
names of the locations within the index and not less than one column 
named 'latitude' and one column named 'longitude'.
Returns
-------
self : object
Instance of the optimizer.
Notes
-----
The gap 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 doesn't have the expected columns"""
if not ('latitude' in df_sites and 'longitude' in df_sites):
raise ValueError("dataframe will need to have columns 'latitude' and 'longitude'")

And now we now have achieved our goal: find the optimal tour from just the sites locations (and never 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 purpose where the optimizer may be very intuitive to make use of. For mere convenience, I’ll add one other method that can be quite helpful in a while after we do [sensitivity analysis] and compare the outcomes of various models. The optimizer, because it is now, tells me the optimal visit order in an inventory, or in a separate dataframe returned by get_tour_stops_dataframe(), but I’d prefer it to inform me the visit order by transforming the locations dataframe that I give it directly—by returning the identical dataframe with a recent column having the optimal sequence of stops. The strategy fit_prescribe can be in command 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 a single line, absorb a dataframe of locations and return 
a replica of it with a recent column specifying the optimal visit order
that minimizes total distance.
Parameters
----------
df_sites : pd.DataFrame
Dataframe with the sites within 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 recent 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 occur, but still
raise Exception("An answer couldn't 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 will solve any TSP in only one line:

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

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

tsp.fit_prescribe(df_sites, sort=False)

And if we’re curious we may also check the variety of variables and constraints the inner model needed to unravel our particular instance of the TSP. This can be handy when doing debugging or performance evaluation.

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