Don Knuth’s MIP, 64 years later

In 1960, the famed computer scientist Don Knuth wrote a technical paper in which he considered an integer programming model for minimizing memory access latency of an IBM 650. I’ve included a screenshot of the 51 variable, 43 constraint problem at the end of this post.

Prof Knuth describes the problem and his attempts to solve it in an unpublished note on his website. Knuth used a technique of Ralph Gomory published the same year, which he somehow became aware of despite the absence of email, the web, arxiv, or large language models. I asked ChatGPT to render an image of Gomory and Knuth pensively sitting atop a mainframe. It refused for privacy reasons. Here is the best it could do:

Knuth ran Gomory’s algorithm on an IBM 650 with less than 10K of memory, but was unable to find an optimal solution. Thirty-five years later (when I was an undergrad), Dimitris Alevras successfully found the optimum value of 22996 using CPLEX on a SPARCstation 5. I am proud to say I have coded on such a machine! Alveras reported that “32,200 branch-and-bound nodes” were needed to solve the problem to optimality. No CPU time is reported in Knuth’s writeup, but we can safely assume the solution time was minutes or likely hours.

I grabbed the problem formulation from Knuth’s writeup and translated it into Python code (copilot was helpful for this purpose). An MPS file with the model can be found here.

I found that the open source solver SCIP was able to find an optimal solution in two tenths of a second, and that Gurobi was able to solve it in less than 0.1s.

This is another testament to the amazingly wonderful effectiveness of operations research!

Finding Optimal State Capitol Tours on the Cloud with NEOS

My last article showed you how to find an optimal tour of all 48 continental US state capitols using operations research. I used the Python API of the popular Gurobi solver to create and solve a traveling salesman problem (TSP) model in a few seconds.

In this post I want to show you how to use Concorde, the world’s best TSP solver for free on the cloud using the NEOS optimization service. In less than 100 lines of Python code, you can find the best tour. Here it is:

Using NEOS is pretty easy. You need to do three things to solve an optimization problem:

1. Create a NEOS account.
2. Create an input file for the problem you want to solve.
3. Give the input file to NEOS, either through their web interface, or by calling an API.

Let’s walk through those steps for the state capitol problem. If you just want to skip to the punchline, here is my code.

Concorde requires a problem specification in the TSPLIB format. This is a text based format where we specify the distances between all pairs of cities. Recall that Randy Olson found the distances between all state capitols using the Google Maps API in this post. Here is a file with this information. Using the distances, I created a TSPLIB input file with the distance matrix – here it is.

The next step is to submit the file to NEOS. Using the xmlrpc Python module, I wrote a simple wrapper to submit TSPLIB files to NEOS. The NEOS submission is an XML file that wraps the contents of the TSPLIB data, and also tells NEOS that we want to use the Concorde solver. The XML file is given to NEOS via an XML-RPC call. NEOS returns the results as a string – the end of the string contains the optimal tour. Here is the body of the primary Python function that carries out these steps:

def solve_tsp_neos_concorde(dist):
    xml = make_neos_concorde(dist)
    neos = NeosClient()
    result = neos.run(xml)
    return tour_from_neos_concorde_result(result)

When I run this code, I obtain the same tour as in my initial post. Hooray! You can also extend my code (which is based on NEOS documentation) to solve many other kinds of optimization models.

Computing Optimal Road Trips Using Operations Research

Randy Olson recently wrote about an approach for finding a road trip that visits all 48 continental US state capitols. Randy’s approach involves genetic algorithms and is accompanied by some very effective visualizations. Further, he examines how the length of these road trips varies as the number of states visited increases. While the trips shown in Randy’s post are very good, they aren’t quite optimal. In this post I’d like to show how you can find the shortest possible road trips in Python using the magic science of operations research. I suggest you read Randy’s post first to get up to speed!

An “optimal road trip” is an ordering of the 48 state capitols that results in the smallest possible driving distance as determined by Google Maps. This is an example of what is known as the Traveling Salesman Problem (TSP). In his work, Randy has made a couple of simplifying assumptions that I will also follow:

• The driving distance from capitol A to capitol B is assumed to be the same as from B to A. We know this isn’t 100% true because of the way roads work. But close enough.
• We’re optimizing driving distance, not driving time. We could easily optimize “average” driving time using data provided by Google. Optimizing expected driving time given a specified road trip start date and time is actually pretty complicated given that we don’t know what the future will bring: traffic jams, road closures, storms, and so on..

These aren’t “bugs”, just simplifying assumptions. Randy used the Google Maps API to get driving distances between state capitols – here’s the data file. Note that Google Maps returns distances in kilometers so you’ll need to convert to miles if that’s your preference.

Randy’s approach to solve this problem was to use a genetic algorithm. Roughly speaking, a genetic algorithm starts with a whole bunch of randomly generated tours, computes their total distances, and repeatedly combines and modifies them to find better solutions. Following the analogy to genetics, tours with smaller total distances are more likely to be matched up with other fit tours to make brand new baby tours. As Randy showed in his post, within 20 minutes his genetic algorithm is able to produce a 48 state tour with a total length of 13,310 miles.

It turns out that we can do better. An inaccuracy in Randy’s otherwise great article is the claim that it’s impossible to find optimal tours for problems like these. You don’t have to look at all possible 48 city road trips to find the best one – read this post by Michael Trick. What we can do instead is rely on the insanely effective field of operations research and its body of 50+ years of work. In an operations research approach, we build a model for our problem based on the decisions we want to make, the overall objective we have in mind, and restrictions and constraints on what constitutes a solution. This model is then fed to operations research software (optimization solvers) that use highly tuned algorithms to find provably optimal solutions. The algorithms implemented solvers rule out vast swaths of possible tours in a brutally efficient manner, making the seemingly impossible routine.

The best-in-class TSP solver is Concorde, and is based on an operations research approach. You don’t need Concorde to solve this TSP – a 48 city road trip is puny by operations research standards. I have chosen to use the Gurobi solver because it is a very powerful solver that includes an easy-to-use Python API, and it has a cloud version. Gurobi even includes an example that covers this very problem! The key to their model is to define a yes-no decision variable for each pair of state capitols. A “yes” value for a decision variable indicates that pair of cities is on the optimal tour. The model also needs to specify the rules for what it means to be an optimal tour:

• The shorter the total distance of the tour (which is determined by the distances between all of the “yes” pairs of cities), the better. This is the objective (or goal) that we seek to optimize.
• The traveller will arrive at each capitol from another capitol, and will leave for another capitol. In other words, exactly two decision variables involving a capitol will be “yes”.
• Tours don’t have cycles: visiting Boise more than once is not allowed.

(Sidebar: If you are not used to building optimization models then the code probably won’t make much sense and you may have no idea what the hell the Gurobi tutorial is talking about. No offense, Gurobi, the tutorial is very well written! The challenge of writing optimization models, which involves writing down precise mathematical relationships involving the decision variables you want solved, is what prevents many computer scientists and data scientists from using the fruits of operations research more often. This is especially the case when the models for classical problems such as the TSP require things like “lazy constraints” that even someone experienced with operations research may not be familiar with. I wrote about this in more detail here. On the other hand, there are a lot of great resources and tutorials out there and it’s simply good practice to rely on proven approaches that provide efficient, accurate results. This is what good engineers do. Anyway, the point is that you can easily steal Gurobi’s sample for this problem and replace their “points” variable with the distances from the data file above. If I wanted to do this with an open source package, or with Concorde itself I could have done it that way too.)

My code, based on Gurobi’s example, is able to find a tour with a total length of 12930 miles, about 380 miles shorter than the original tour. What’s more, it takes seconds to find the answer. Here is my Python code. Here is the tour – click here to explore it interactively.

A text file with the tour is here and a GPX file of the tour is here courtesy of gpsvisualizer.com. This optimal tour is very close to the one the genetic algorithm came up with. Here is a screenshot for reference:

An interesting twist is that Randy extends the problem to consider both the driving distance and the number of states visited. If we are willing to do a tour of, say, 10 states, then clearly the total distance for the tour will be much shorter than a 48 state tour. Randy has a nice animation showing tours of differing numbers of states, as well as a chart that plots the number of states visited against the estimated driving time. This curve is called the efficient frontier – you sometimes see similar curves in financial models.

The modified problem of finding the shortest tour involving K of 48 total state capitols can also be solved by Gurobi. I extended the optimization model slightly:

• Introduce new yes-no decision variables for each capitol: “yes” if the capitol is one of the lucky K to be visited.
• Exactly K of the new decision variables should be “yes”.
• Fix up the original model to make sure we don’t worry about the other N-K cities not on our mini tour.

(I also had to modify the model because I am running on the cloud and the “lazy constraints” mentioned in Gurobi’s TSP writeup don’t work in the cloud version of Gurobi.)

With this new code in place I can call it for K=3…47 and get this optimal efficient frontier curve:

The distances and tours for all of these mini tours are given here.

What have we learned here? In about 200 lines of Python code we were able to efficiently find provably optimal solutions for the original road trip problem, as well as the “pareto optimization” extension. If you’re a data scientist, get familiar with operations research principles because it will certainly pay off!

Detecting the Sources of Model Infeasibility using Gurobi

Today I come to sing the praises of Irreducible Infeasible Sets (IIS). An IIS is a minimal subset of the constraints of an optimization problem that are self-contradictory. In other words, if you remove any single constraint in an IIS, then the remaining constraints are feasible. Most modern linear and mixed-integer solvers have programming APIs to return IIS information for an infeasible model.

The reason IIS is so cool is that it helps diagnose why a model is not able to be solved. This is especially useful when the optimization problem to be solved is formed with the help of input from a user who may not know anything about operations research, or when the model is complicated. IIS is like the “spell check” of linear programming. Building optimization models without IIS is like debugging using “printf”.

Here is a simple example from the lp_solve documentation:

min x + y
  x >= 6
  y >= 6
  x + y <= 11

6 plus 6 is more than 11, people. In this case, all three constraints are part of the IIS.

Here is some C# code using the Gurobi API that obtains the IIS for a very simple model. First, here is the code that forms and solves the model. This is easy to figure this out by reading the Gurobi docs. The only twist is that I have expressed the first two constraints as lower bounds on the variables x and y.

GRBEnv env = new GRBEnv();
GRBModel model = new GRBModel(env);
GRBVar x = model.AddVar(6, Double.MaxValue, 1, GRB.CONTINUOUS, "x");
GRBVar y = model.AddVar(6, Double.MaxValue, 1, GRB.CONTINUOUS, "y");
model.Update();
GRBConstr c3 = model.AddConstr(x + y, GRB.LESS_EQUAL, 11, "c3");
model.Optimize();
model.ComputeIIS();

To find the IIS you need to call another method called ComputeIIS. Then you can query the constraints and variables to see if they are part of the IIS. This is accomplished using the Get method. I simply print out the constraint and variable names that are in the IIS.

// Print the names of all of the constraints in the IIS set.
foreach (var c in model.GetConstrs())
{
    if (c.Get(GRB.IntAttr.IISConstr) > 0)
    {
        Console.WriteLine(c.Get(GRB.StringAttr.ConstrName));
    }
}                

// Print the names of all of the variables in the IIS set.
foreach (var v in model.GetVars())
{
    if (v.Get(GRB.IntAttr.IISLB) > 0 || v.Get(GRB.IntAttr.IISUB) > 0)
    {
        Console.WriteLine(v.Get(GRB.StringAttr.VarName));
    }
}

If you want to write the IIS information to a text file for review, simply use the write method with the “ilp” extension on the file name:

model.Write("model.ilp");

One last thing: using LINQ (a feature of the C# language), you can simplify how we query for the IIS constraints and variables. For example, I can get all of the constraint names in one statement as follows:

var allConstraintNames =
    from c in model.GetConstrs()
    where c.Get(GRB.IntAttr.IISConstr) > 0
    select c.Get(GRB.StringAttr.ConstrName);

LINQ is very handy for creating models and querying results. I use it often!