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’ve got the simpler 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 wish the optimizer to seem like a “magic box” that even users not acquainted with 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:
- 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 vital difference between what we did and what we’re doing is that, now, the initial site just isn’t hard-coded as simply"hotel"
, but is assumed to be the positioning of the primary row indf_distances
. In the overall 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. - (Try and) Solve the model. That is performed within the
_optimize
method inherited fromBaseOptimizer
, which returnsTrue
provided that an answer is found. - Extract the answer from the model and parse it in a way that is straightforward 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 attributestour_
andtour_distance_
, respectively. This method gets invoked provided that an answer exists, so if no solution is found, the “solution attributes”tour_
andtour_distance_
never get created. The good thing about 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 might be conveniently retrieved at any point, not necessarily in the intervening time 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 may 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 wish 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 incorporates every time we print it.
from typing import Tuple, Listclass TravelingSalesmanOptimizer(BaseOptimizer):
"""Implements the Miller–Tucker–Zemlin formulation [1] of the
Traveling Salesman Problem (TSP) as a linear integer program.
The TSP might be stated like: "Given a set of locations (and typically
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 provide 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 just isn't 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 have 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 after all 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. For those who 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’d have the reply to our problem. We don’t must 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 will be enough to make use of the rank variables to extract the optimal tour. Nevertheless, because the TSP is just the initial prototype of what is going to grow to be 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 following 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, regardless of how complex the model we’re using.
The advantages of this approach will grow to be apparent in the next articles, where latest 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 have now the variable 𝛿ᵢⱼ, indexed by possible arcs (i, j), where 𝛿ᵢⱼ=0 means the arc just isn’t chosen and 𝛿ᵢⱼ=1 means the arc is chosen.
- We would like 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):
- Extract the chosen arcs from 𝛿ᵢⱼ, which represents the tour. Done in
_get_selected_arcs_from_model
. - Convert the list of arcs (edges) into a listing of stops (nodes). Done in
_get_stops_order_list
. - 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 a listing of ordered stops just isn’t that straight-forward. To avoid convoluted code, we exploit the undeniable 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 a listing **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 latest 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’ve got the answer in a format easy to interpret and work with. However the sequence of stops just isn’t all we care about. The worth of the target function can also be a very important metric to maintain track of, as it is the criterion guiding our decisions. For our particular case of the TSP, this implies getting the whole 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 whole distance of the tour, regardless that, here too, each measures are equivalent. In additional complex models, the target function will even incorporate other targets, so this equivalence will now 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 small print 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 whole distance of the optimal tour"""
if self.is_fitted:
# as the target is an expression for the whole distance,
distance_tour = self.model.objective() # we just get its value
return distance_tour
print("Optimizer just isn't fitted to any data, no optimal objective exists.")
return None
It might look superfluous now, nevertheless 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 vital, 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 latest 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 can 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 will 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] incorporates 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 are able to 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 give you 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 no less than one column
named 'latitude' and one column named 'longitude'.
Returns
-------
self : object
Instance of the optimizer.
Notes
-----
The space 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 should have columns 'latitude' and 'longitude'")
And now we’ve got 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 could 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 a listing, 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 brand new column having the optimal sequence of stops. The tactic fit_prescribe
can be accountable for 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, soak up a dataframe of locations and return
a replica of it with a brand new 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 brand 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 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 are able to solve any TSP in only one line:
tsp = TravelingSalesmanOptimizer("Paris")tsp.fit_prescribe(df_sites)
If we’d wish to conserve the unique order of locations as they were in df_sites
, we are able to 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 interior 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