We need better tools for linear programming—here's how and why
The Berlin airlift, GAMS, AMPL, Pyomo, JuMP, Linopy, gurobi-pandas and more
In 1949, at the start of the cold war, the Soviet Union blocked all land access to Berlin in a move to pressure Western Allies. Unwilling to seed to Soviet demands, American and British forces undertook a humongous logistical operation to sidestep the blockade. For over a year, they sent ~850 flights per day into West Berlin to deliver food, coal and supplies in what is known as the Berlin Airlift.
Back home in the United States, scientist George Dantzig was trying to develop a process to reliably find solutions to the Pentagon’s logistics problems. He discovered that if he expressed a logistical problem as a set of linear equations he could then use a mathematical algorithm he had invented to find the optimal logistical plan1. He used the Berlin Airlift logistics to demonstrate the usefulness of his algorithm and in doing so solved the first ever linear program. The field of convex optimization was born! [1]
Today, dozens of sectors rely on convex optimization. Which stocks should a hedge fund buy to make a profit-maximizing risk-minimizing portfolio? How should an airline allocate pilots and flights to meet customer demand while minimizing costs? Where should a city build their next subway line to serve the most commuters? Convex optimization problems are everywhere.
I specialize in energy systems engineering where convex optimization is at the core of our quest to building net-zero electricity grids. Using cost and performance estimates for technologies like solar panels, batteries and transmission lines, convex optimization allows us to run large simulations to answer important questions. When might hydrogen be a better energy storage technology than batteries? How much energy storage do we even need? How might nuclear power help us decarbonize electricity grids? Convex optimization is the best tool we have for predicting our inherently unpredictable future.
Despite convex optimization being central to so many fields, I’ve felt deeply unsatisfied with the tools for convex optimization in my work as a researcher and professional power systems analyst. So, last semester I recruited my friend Emily Traynor and together we set off to make things better. Since then, I’ve spoken to several professors, researchers, and industry professionals to understand the state of convex optimization tools. Here’s what I’ve learnt.
Great solvers, poor solver interfaces
In convex optimization there are two types of tools: solvers and solver interfaces.
Solvers, such as Gurobi, are highly optimized pieces of software that solve convex optimization problems using mathematical algorithms. Solver performance matters a lot since problems, often with tens of millions of variables, can take hours or days to solve. Despite no major algorithmic breakthroughs since the mid-1980s, the performance of solvers is continuously improving due to increases in computing power and the development of highly effective heuristics and code optimization tricks. For example, a recent release of the Gurobi solver boasted a 25% speedup and the capability to solve previously intractable problems.
Gurobi’s focus on performance, as well as its free-for-academia policy, has led it to dominate the solver market. Typically, only non-academic users who aren’t willing to pay $10 to $50k / year for Gurobi use other solvers like CPLEX or the open-source HiGHS. I don’t expect Gurobi’s dominance to change anytime soon unless a competitor makes an algorithmic breakthrough (e.g. parallelization of solving on GPUs).
Unfortunately, a solver alone is not of much use since solvers operate on optimization problems in their mathematical form—large matrices that look like gibberish to humans.
To generate these matrices practitioners need to use a solver interface. A good solver interface allows practitioners to define an optimization problem at a high level—at the level of the objective, constraints and decision variables of their real-world problem—while handling all low-level details like generating the matrix and sending it to a solver.
Compared to solvers which have seen continuous improvement and fierce competition in recent years, solver interfaces have seen little progress.
The fall of algebraic modelling languages
In 1978, GAMS was created and became the first modern algebraic modelling language. Like programming languages, GAMS’ users write code to define optimization problems, but unlike programming languages GAMS code defines algebraic equations, not a sequence of executable steps.
In 1985, AMPL, another algebraic modelling language, was created to improve on GAMS’ non-intuitive syntax [2]. Today, GAMS and AMPL are the predominant solver interfaces used in industry despite costing ~$7k/year. Their success stems from the inherent properties of algebraic modelling languages—high performance and a syntax close to mathematical notation that makes for terse and clear optimization models.
Unfortunately, algebraic modelling languages like GAMS and AMPL are failing to keep up with the modern software ecosystems. Neither language integrates well into modern IDEs, supports autocomplete, or has a code formatter. To make matters worse, GAMS uses its own propriety .gdx datafile format. One policy analyst I spoke to described GAMS as “archaic.”
This archaicness, costs and learning curve have caused algebraic modelling languages to fall out of favour in light of newer open-source solver interfaces, despite the inherent performance and syntax benefits of algebraic modelling languages.
The challenges with open-source solver interfaces
In 2008, researchers at Sandia National Laboratories working on “national security applications” developed Pyomo2, an open-source solver interface built entirely in Python. Instead of needing a custom algebraic modelling language, Pyomo users would define their optimization models directly within their Python code which allowed for much more flexibility.
Although flexible and free, Pyomo is >10x slower at constructing large models compared to GAMS or AMPL [3]. So, in 2023, frustrated researchers at TU Berlin published Linopy3, a Python solver interface that is much faster than Pyomo due to its use of vectorization4. Unfortunately, Linopy doesn’t work with very sparse datasets that are common in many optimization problems.
Another notable attempt at building a great open-source solver interface came out of MIT in 2017. Researchers built JuMP in the programming language Julia. JuMP has a particularly clear syntax due to its use of syntactic macros, a feature of Julia that doesn’t exist in Python5. Unfortunately, most users aren’t willing to switch from Python to Julia.
Gurobi, the undesired solution
The solver interfaces discussed up to now have all been solver agnostic. For example, convex optimization models defined in Pyomo can be solved using Gurobi, CPLEX, HiGHS or any one of the many other existing solvers.
However, there exists solver-specific solver interfaces. For example, Gurobi offers a Python API that wraps its solver. Users can call this API to construct their models. In fact, in 2022 Gurobi released the library gurobi-pandas to make creating models with Gurobi particularly easy. Compared to my work with GAMS, Pyomo and Linopy, gurobi-pandas provided the best experience and performance so far.
Unfortunately, defining your model in a Gurobi specific solver interface forces you to use Gurobi indefinitly. This is particularly problematic for open-source models used by multiple groups since not everyone has access to the Gurobi solver.
The future of solver interfaces
Every solver interface to date fails at least one important requirement. So what will the next solver interface look like?
Most likely, the next solver interface will be a Python library using vectorization like Linopy. However, instead of relying on multidimensional arrays that make it inherently difficult to work with sparse datasets, the next solver interface could use simple 1D arrays via Pandas or Polars.
Other options exist too. AMPL’s main failure is being costly and closed source. Perhaps the syntax and performance benefits of algebraic modelling languages could be revived if an open-source AMPL-like algebraic modelling language was built. Unfortunately, this is no easy task for the open-source community.
Finally, when my friend Emily and I decided to explore this space, we were hoping to build a startup. Perhaps the next solver interface will allow users to build optimization models by simply typing out their mathematical equations, without the need to write any code. It wouldn’t be the first time startups have tried to build no-code solutions.
References
[1] An interview with George Dantzig, father of linear programming
[2] An interview with Brian Kernighan, co-founder of AMPL
[3] JuMP paper
A logistical plan at the time was called a logistical program. This is where the term linear programming comes from.
Why do good things always come for the military :(
Linopy spun out of a project called ‘nomopyomo’. Goes to show how frustrated researchers were with Pyomo.
Vectorization is done using the xarray library.
There’s a proposal to add syntactic macros to Python but it hasn’t gone anywhere.