Everything started in the middle of a rainy Saturday, I was binge-watching one of the new Christmas movies on Netflix when I heard someone loudly shouting my name just outside my apartment. I opened the door and I saw a big figure wearing a red coat, red pants, and a funny red hat. I’ve never been so surprised in my life: the man staring at me was Santa in person!!!

TL;DR; what he wanted from me was a piece of software to optimize his delivery for December 24th. His request was to find the shortest route given 342 cities. He needed a cyclic path starting and ending in Lapland.

Here’s a video for convenience! The link to the live JavaScript visual is at the end of the article!

What follows is the conversation I had with Santa!

The Santa’s Problem

Andrea: Hi Santa, let me prepare a cup of tea, in the meanwhile, why don’t you explain to me your problem?

Santa: Ho Ho Ho! Well, as you know every year during the night between the 24th and 25th of December, I need to deliver billions of toys all over the world. It’s a very challenging task, Andrea. I’m wondering if Computer Science can simplify my life. How can visit every city travelling for the shortest distance possible?

Andrea: Ok, this sounds familiar, can you expand a bit on your problem?

Santa: Of course. I have a list of 307 cities, those with the largest population. I would like to know which is the shortest path that visits all of these cities only once. As you might know, I need to depart from my house in Lapland, and Lapland is also the last stop of my trip! Obviously, I need to touch more than 307 cities, but I want a strategic plan, not operational. Once I’ve reached the major places, I’ll figure out with the reindeers how to move to the nearby towns!

Andrea: Ok Santa. I think your request is clear. Actually what you are describing is a very classical problem in Computer Science: it is called the Travelling Salesman Problem (TSP), and its qualitative formulation is pretty much what you just told me. I have bad news though: it’s very difficult to solve! Let me draw a small sketch to see if we agree on what’s the task!

Let’s say that we want to focus only on the major US cities, then the route that you are looking for is something like the following:

Santa: Oh Andrea that’s exactly what I am asking for! Just remember that I can fly! By the way, can you explain to me why it is so difficult?

Andrea: Well, the TSP belongs to a class of problems known as NP-Complete problems. Not sure you want to know the mathematical definition, but from a practical standpoint, in over a hundred years of research, nobody has found a fast method to solve any of those and these solutions are unlikely to exist.

Santa: How hard can it be?

Andrea: If you run a brute force algorithm for what you are asking (which means testing all the possible solutions and picking the best one), it can take LITERALLY hundreds or thousands of years to get the results! Even with modern technologies. Luckily for you, there are some smarter ways to deal with many TSP instances!

Santa: Ho Ho Ho! This is great news! I am happier than that day I fell in the marshmallows pool! Can you tell me what are the alternatives?

Andrea: Of course. There are a few strategies:

If you just want a “good enough” solution, we could use heuristic algorithms : these methods are usually really fast, but they lead to solutions that may not be (and usually are not) optimal .

: these methods are usually . We can use some smart “exact” algorithm : potentially it will take forever, usually they solve small problems in an acceptable time , but sometimes they are difficult to implement. In this case, the solution is always optimal.

: potentially it will take forever, usually , but sometimes they are difficult to implement. In this case, the solution is always optimal. We can use Mixed-Integer Linear Programming: if we follow this route, we only need to write a mathematical model and then use a MILP solver as a black box to get the optimal solution.

Before recommending you something, may I know what is your budget and what is your deadline?

Santa: Oh well, unfortunately, my deadline is quite strict, I need the solution by Sunday night, so you would have a day and a half to implement everything. On the other hand, I have an infinite budget! The National Bank of Christmas can print as many Magic Dollars as I want! Note also that I want nothing else than the optimal route among the 307 cities!

Andrea: Ehm… ok. Indeed the deadline is quite close; since you need the optimal route we are going to go with Mixed-Integer Linear Programming (MILP). There are many tools that we can use, some of them are commercial, like CPLEX or Gurobi. Some of them are free to use only for academics, like SCIP. Some of them are open-source, such as GLPK and CBC.

Now, you told me that you have an infinite budget, but unfortunately, Gurobi and CPLEX won’t accept Magic Dollars as currency, so we’ll need to use GLPK or CBC; the latter is usually faster and can be used with Python, so I suppose we have our base tools!

Santa Modeling

Santa: Wooow! Amazing Andrea! Let’s start then. Many are unaware that a long time ago I graduated in Mathematics. I would be really happy if you’d explain the model.

Andrea: Oh great, it’s always fun to write down some formulas, it makes me feel smarter! First, let’s see the solution that we are looking for in a sample situation:

This can be roughly described as follows:

It is optimal , indeed if we consider all the possible connections among the nodes, we only selected those that drive to the shortest cycle that visits every point.

, indeed if we consider all the possible connections among the nodes, we only selected those that drive to the shortest cycle that visits every point. Each city is visited only once , this means that we enter only once in each city and we exit once from each city ( more formally we could say that each node has exactly one incoming edge and exactly one outgoing edge ).

, this means that we enter only once in each city and we exit once from each city ( ). We only have one single cycle!

To express this solution using Linear Programming we need to define some variables: given that we have n cities, we could create one variable for each pair of cities i and j, yᵢⱼ.

Now, yᵢⱼ will be equal to 1 if in our best route we need to go from city i to city j straight forward (directly, without first passing through another city), yᵢⱼ will be equal to zero in case that connection is not chosen. cᵢⱼ is the cost of the connection between city i and city j. Let me write the model:

Santa, the issue with this object is about the “No Sub-tours” constraints (e.d. see the image above). A sub-tour is a cycle that does not visit all the cities, but only a subset. We want to have a single cycle, so those inequalities aim to remove the possibility of having more than one cycle.

To implement this model we should write roughly one “No sub-tour” constraint for each subset of nodes. If we have n nodes that would mean 2ⁿ-2 constraints which are impossible to represent in the computer’s memory (e.g. with n=300, the number of constraints would be 20.370.359.763.344.860.862.684.456.884.093.781.610.514.683.936.659.362.506.361.404.493.543.812.997.633.367.061.833.97.376 minus 2).

This is why we need to find a better strategy.

What we can do is apply the following algorithm:

We remove the “No Sub-tours” constraints from the model, obtaining what is called a “relaxed model”. We solve the relaxed model to the optimum. This solution will very likely contain sub-tours. We find all the sub-tours in the current solution and we add to the model those constraints that forbid those sub-tours. We go to point 2 and iterate until no sub-tours are detected.

Maybe later we should make a sketch to clarify these steps.

Santa: Oh, that’s fantastic Andrea! Thanks for the explanation! So the idea is to remove that huge number of constraints and add only those that are necessary to obtain a solution with a single cycle. Ahhh sweet memories! And is this the best method to solve this problem?

Andrea: Of course not Santa. There are many other types of constraint that we could add to the model to reduce the number of iterations needed to get to the optimal solution (they are called “cuts”, from cutting). But you told me that we only have two days so we’re going for the simplest way.

Drawing the algorithm

Santa: Ho Ho Andrea, I would also like to see some code. You know, the whole IT System that we have at Santa’s House is written in Pascal, we’d really need some renovation. Could we do some old-but-gold pair programming so that I can start to delve into Python?

Andrea: Oh wow, absolutely Santa. Let’s implement this together.

As I told you, we need to relax the model removing “No sub-tour” inequalities (these are technically called SECs, Sub-tour Elimination Constraints). Our model becomes the following:

As I explained earlier, if we solve the above we may have those very famous sub-tours. For example, the following image is the solution to the relaxed model using 18 cities:

As you can see the solution is satisfying all the constraints: it is optimal and each city has only one incoming edge and one outgoing edge. But we do not want sub-cycles! As I said, we need to:

Identify all the “Connected Components” (intuitively all the groups of nodes that are connected together). Update the model with all the inequalities that forbid those components.

Let me draw a diagram:

Now that we have added all the new cuts, let me show you what happens if we solve the problem:

We begin to see the shape of the solution! Finally, if we iterate the process… Boom!

The Implementation

Andrea: It’s time to write some code. First of all, we are going to need some libraries. The most important one is PuLP. From their website: PuLP is an LP modeler written in Python. PuLP can generate MPS or LP files and call GLPK, COIN CLP/CBC, CPLEX, and GUROBI to solve linear problems. PuLP comes bundled with CBC. We also need networkx, a Python package for the creation, manipulation, and study of the structure of networks (we need it to identify the sub-tours during the optimization process).

Let’s create a Python class that will solve the problem. First, we need to read the cities (here we can have the list of the world’s cities with the population and the geographical coordinates). Since your house location is not in the list, Santa, we’ll add it to the data; moreover, we will add some more cities to improve the world’s coverage:

Now, we need a function to calculate the distance between two locations (using latitude and longitude coordinates). Your vehicle is a sleigh pulled by magical reindeers, so the distance as the crow flies will be a good approximation:

Nice! We can write the function to build the model:

Good. It’s time for the complicated part.

We have to get the current solution and build the corresponding graph. We’ll use networkx. After that, we need a function to:

Detect the connected components

If those are more than 1 (meaning that we have sub-tours), add the new constraints to the model. If we only have one connected component, then we’ve found the optimum.

We are ready to write the function that solves the model!

Done, dear Santa, we just have to call the methods!

Santa: I can’t wait to execute the code!

Andrea: Me too! Let’s see how long it takes and how many SEC constraints we need to add to get to the solution! We’ll use JavaScript to visualize the results; more specifically we’ll implement the visuals on a map using Leaflet along with Leaflet.motion plugin!

I’ll run the code on my laptop:

MacBook Pro (13-inch, 2019)

Processor: 2.8 GHz Quad-Core Intel Core i7

RAM 16 GB 2133 MHz LPDDR3

First of all, we will try to use your 308 cities (307 plus your house in Lapland):