Some notes on the landscape of linear optimization software and applications

Kragen Javier Sitaker, 2019-08-21 (updated 2019-08-25) (35 minutes)

I started auditing an “operations research” class at the UBA this quatrimester. It’s largely focused on linear programming, which is to say, optimizing linear functions subject to linear constraints, but the class also covers some other non-linear-programming topics.

I’m particularly interested in optimization algorithms, as I’ve said a few times before, because they potentially raise the level of abstraction in programming — rather than specifying how to compute something, you just specify how to recognize if something is the thing you wanted to compute, and then let some kind of generic solver figure out how to solve the problem. That’s what earned them such a prominent place in More thoughts on powerful primitives for simplified computer systems architecture. And it turns out that linear-optimization solvers, as they’re called, are the most developed family of optimization algorithms out there.

Other useful overviews of the space can be found in the NEOS Server Optimization Guide and the GLPK Wikibook.

I’m paying as little attention as possible to proprietary software here, but unavoidably need to mention it periodically simply because it has played such a central role in this space.

Overview of linear programming

Linear programming is a particularly popular subset of mathematical optimization because there are well-known algorithms for solving it that scale to fairly large numbers of decision variables (the variables the optimization algorithm is allowed to decide on a value for to reach the optimum), and it provides a useful approximation for many economically important problems. It’s so popular that sometimes it overshadows the rest of mathematical optimization — many books purportedly about mathematical optimization consist almost entirely of material on linear optimization, with perhaps a short section at the end on nonlinear optimization.

I’ve tended to be a bit prejudiced against linear optimization because the real world is never perfectly linear and because I felt that linear constraints weren’t very expressive. But, now that I’m taking this class, I’m starting to appreciate the astounding expressiveness of linear programming — especially “mixed integer” linear programming (“MIP” or occasionally “MILP”), an extension of pure linear programming in which some of your decision variables are constrained to integer values, which is more expressive but enormously more difficult to solve.

As the lp_solve documentation summarizes it:

The linear programming (LP) problem can be formulated as: Solve A·x >= V₁, with V₂·x maximal. A is a matrix, x is a vector of (nonnegative) variables, V₁ is a vector called the right hand side, and V₂ is a vector specifying the objective function.

An integer linear programming (ILP) problem is an LP with the constraint that all the variables are integers. In a mixed integer linear programming (MILP) problem, some of the variables are integer and others are real.

(Some formulations additionally specify that all xᵢ are nonnegative, which turns out to matter very little.)

There are several algorithms known for solving linear optimization problems, of which the oldest and most-widely-used is George Dantzig’s “simplex algorithm”, despite its worst-case exponential complexity. Many newer solvers instead use the “interior-point method”, which has worst-case polynomial complexity for continuous linear optimization. (Most mature solvers support both.) However, for mixed integer linear programming, they all must fall back on exponential-time search procedures, usually “branch-and-cut” procedures (also known as “branch-and-bound”).

The software landscape

Typically, as an end-user, you solve linear optimization problems by writing your problem in an “algebraic modeling language” such as GNU MathProg or “GMPL” (a dialect of a popular proprietary modeling language called “AMPL”) or ZIMPL; the model translator compiles your model (including data files, which it might pull from a database) into a standard format (such as “MPS”, “LP”, or .nl — see below), generally expanding it considerably, and submits it to a “solver”.

Unfortunately, the most popular model translators and solvers are all proprietary; proprietary CPLEX, now owned by IBM, is the most popular solver, with proprietary Gurobi and proprietary SCIP second and third, and proprietary AMPL is the most popular model translator; there is an excellent AMPL book by Brian Kernighan (et al., but “Kernighan” is a sufficient recommendation) which also serves as an introduction to linear programming. Both CPLEX and SCIP are gratis for institutional academic use (we are using SCIP in the class) but not to independent researchers; I don’t know about AMPL. So we can forget about these.

On the NEOS server, anyone can use most of this software in batch mode for free, for optimization jobs of up to 16 megabytes of input, 3 GB of RAM, and 8 hours of run time, thanks to the Wisconsin Institute for Discovery at UW Madison; the NEOS server has AMPL, ZIMPL, CPLEX, SCIP, CLP, CBC, Gurobi, CBC, and SYMPHONY, but not including GLPK or lp_solve, perhaps because it doesn’t consider them “state of the art”.

XXX I really need a diagram!

Modeler–solver interfaces (MPS, Osi, .nl, and LP)

Often the model translator talks to the solver through a file in some standardized file format; MPS is the most popular such format, a punched-card-based format originally designed for IBM’s Mathematical Programming System/360. LP, sometimes called “CPLEX LP” because of its origins in CPLEX, is a more-readable but less-widely-supported format; and AMPL, a popular but proprietary algebraic modeling language (compatible with GLPK’s MathProg), defined its own .nl format, mentioned earlier, that has been interfaced to a wide variety of solvers.

Linear-programming solvers have been around considerably longer than algebraic modeling languages, so both MPS and LP are designed for humans to write, but MPS is designed badly.

(There are actually two different, incompatible MPS formats: the original fixed-field MPS format from the 1960s, and a slightly less bletcherous “free-form” 1980s version.)

Here’s a small LP file that GLPK and CBC are able to accept and solve (the optimum is z = 1, y = 0.33...).

Maximize
 x: + 2 y + 3 z
Subject To
 a: + 2 z + 3 y <= 3
Bounds
 0 <= y <= 4
 0 <= z <= 1
End

You can save that in foo.lp and solve it with GLPK by running glpsol --lp foo.lp -o foo.out. You can transate it to MPS with glpsol --check --lp foo.lp --wmps foo.mps.

Most solvers, including GLPK’s, also have bindings that allow you to call them from your program without preparing an input file for them to parse; and the COIN-OR project has written a standard API called “Osi” for doing this with many different solvers.

GLPK (modeling language (GMPL/MathProg) and solver)

GLPK, a C library, is the thousand-pound gorilla of linear optimization systems; Octave is built with it, and there are bindings to it in Java, Python, and R. The GLPK package includes the model translator for the popular algebraic modeling language GNU MathProg (aka GMPL), and also solvers for the revised simplex method, the primal-dual interior point method, and the branch-and-bound method.

MathProg is a subset of the AMPL language mentioned above, including nearly all the facilities of AMPL for defining optimization problems, but without statements like include, model, and data.

It supports MPS and LP format files for both input and output, and can even translate between them, as well as generating them with the MathProg model translator.

GLPK was originally released in 2000, but has been under active development ever since; there is a Wikibook about it, linked above in the introduction.

Here’s a MathProg version of the example LP file above, with additional directives to solve it and display the results:

var y >= 0, <= 4;
var z >= 0, <= 1;
maximize x: 2*y + 3*z;
subject to a: 2*z <= 3 - 3*y;
solve;
display x, y, z;
end;

If you save this as min.mod you can run glpsol -m min.mod to get the solution, or glpsol -m min.mod --wlp min.lp to get a slightly expanded version of the LP file above. Pitfall: if it’s gzipped, like most of the examples that come with GLPK, glpsol will uncompress it to run it, but at least the version I have dies with a spurious “read error”.

This example is somewhat freer-form than the LP file (there’s an inequality with decision variables on both sides, for example), but to really show the advantages of an algebraic modeling language like MathProg, we need a bigger model. For example, here’s an integer-linear-programming sudoku solver I wrote this afternoon:

param N default 9; set S := 1..N; param SW 'square width' default 3;
set G 'givens' default {};
param givenrow {G}; param givencol {G}; param givenchoice {G};

var X { S, S, S } binary;

maximize pleasure: 69;

rows {col in S, choice in S}: sum {row in S} X[row, col, choice] = 1;
mutex {row in S, col in S}: sum {choice in S} X[row, col, choice] = 1;
cols {row in S, choice in S}: sum {col in S} X[row, col, choice] = 1;
squares {bigrow in 0..SW-1, bigcol in 0..SW-1, choice in S}:
    sum {row in 1..SW, col in 1..SW}
        X[bigrow * SW + row, bigcol * 3 + col, choice] = 1;
givens {j in G}: X[givenrow[j], givencol[j], givenchoice[j]] = 1;

glpsol --check --wlp sudoku.lp -m sudoku.mod translates this 12-line MathProg model into a CPLEX LP file of 2000+ lines, beginning as follows:

\* Problem: sudoku *\

Maximize
 pleasure: 0 X(1,1,1)
\* constant term = 69 *\

Subject To
 rows(1,1): + X(1,1,1) + X(2,1,1) + X(3,1,1) + X(4,1,1) + X(5,1,1)
 + X(6,1,1) + X(7,1,1) + X(8,1,1) + X(9,1,1) = 1
 rows(1,2): + X(1,1,2) + X(2,1,2) + X(3,1,2) + X(4,1,2) + X(5,1,2)
 + X(6,1,2) + X(7,1,2) + X(8,1,2) + X(9,1,2) = 1
 rows(1,3): + X(1,1,3) + X(2,1,3) + X(3,1,3) + X(4,1,3) + X(5,1,3)
 + X(6,1,3) + X(7,1,3) + X(8,1,3) + X(9,1,3) = 1

There’s a SciTe-based IDE for GLPK called Gusek I haven’t tried. I’ve often had a lot of trouble debugging my MathProg models, and I wonder if it might help.

ZIMPL (modeling language)

ZIMPL, the Zuse Institute Mathematical Programming Language, is a modeling language from the same academic research institute, the Zuse Institute Berlin, that wrote the solver SCIP, but ZIMPL is free software, while SCIP is not. It’s mostly compatible with GNU MathProg. It started out as Thorsten Koch’s Ph.D. thesis in 2004, ZIB-Report 04-58, “Rapid Mathematical Programming”. (“Mathematical programming” is a synonym for “mathematical optimization”.)

The ZIMPL language appears to be more or less at the same level of abstraction as MathProg, but with arguably nicer syntax. So for example in ZIMPL you might say

minimize cost: 12 * x1
    + sum <a> in A : u[a] * y[a];
subto oneness: sum <a> in A : u[a] == 1;
subto nonnegative: forall <a> in A : y[a] >= 0;

while in AMPL/MathProg you would say

minimize cost: 12 * x1
    + sum {a in A} u[a] * y[a];
subject to Oneness: sum {a in A} u[a] = 1;
subject to Nonnegative {a in A}: y[a] >= 0;

It can produce MPS and LP files, as well as something called “Polynomial IP”. Beware, it generates the name for its output file from the name of the input file, and so it will overwrite any existing file of the corresponding name.

Here’s a translation of the tiny model used as examples earlier into ZIMPL:

var y >= 0;
var z >= 0;
maximize x: 2*y + 3*z;
subto a: 2*z <= 3 - 3*y;
subto ym: y <= 4;
subto zm: z <= 1;

If saved as foo.zpl, you can convert it to LP format with zimpl foo.zpl, writing the output to foo.lp.

Since ZIMPL and MathProg are so similar in their capabilities, but MathProg is bundled with the popular GLPK solver, I don’t know what the advantage of using ZIMPL would be — except that the proprietary solver SCIP has ZIMPL built in, and so it’s more convenient to use SCIP with ZIMPL than with MathProg.

The biggest difference I’ve found so far between ZIMPL and MathProg, actually, is that MathProg lets you specify things to do after the solution is found, at least when you’re using the GLPK solver. This way you can present the solution in a more readable form. My Sudoku solver above, for example, has half a page of code tagged onto the end to display an ASCII-art Sudoku board. To do the same thing with SCIP and ZIMPL, I ended up writing a Python script to match regexps against the SCIP output.

CMPL (modeling language)

CMPL is another algebraic modeling language from the COIN-OR project, GPLv3 licensed, bundled with an IDE called Coliop; it can use the CBC or GLPK solvers, as well as proprietary solvers. It also includes Java and Python interfaces, presumably as a sort of embedded DSL. The project mailing list seems consistently very low traffic since 2011.

I haven’t tried it yet, although it looks like it might address some of my annoyances.

MiniZinc (modeling language)

MiniZinc is a constraint modeling language that is not limited to linear constraints; it can use CBC as a solver, as well as MinisatID, proprietary SCIP, and over a dozen other solvers. I haven’t tried it yet.

libisl (ILP solver)

The documentation says:

isl is a library for manipulating sets and relations of integer points bounded by linear constraints. Supported operations on sets include intersection, union, set difference, emptiness check, convex hull, (integer) affine hull, integer projection, and computing the lexicographic minimum using parametric integer programming. It also includes an ILP solver based on generalized basis reduction.

I have it on my laptop because GCC depends on it, apparently for a newish loop optimization framework internal to GCC called Graphite. It sounds scarily powerful. Unlike the others, it’s specifically intended for integer optimization, with no support for continuous-domain optimization.

Needless to say, it doesn’t interoperate with the others.

COIN CLP/CBC/SYMPHONY (solvers)

The COIN project is an Apache project for making operations-research results reproducible by basing them on open-source software. It includes the solvers CLP, CBC, and SYMPHONY; CBC and SYMPHONY require an underlying linear-programming solver like CLP; it also defined an API called “Osi”, the Open Solver Interface, as an alternative to the file-format interfaces described earlier. I haven’t yet figured out how to use SYMPHONY yet.

These three solvers are the only free solvers I’ve found on the NEOS server.

I’ve been able to solve minimization MPS problems with CLP using clp foo.mps but solving maximization problems or seeing the actual results requires a slightly more elaborate command:

$ clp solar.mps -maximize -dualsimplex -solution solar.solution

If you apply this to an MPS file containing an integer-programming problem like the Sudoku solver given above, you get garbage, because CLP finds a continuous solution instead of an integer solution. There is a similar cbc executable which has an -branch option to avoid this:

$ cbc sudoku.mps -dualsimplex -branch -solution sudoku.solution

This finds an answer which appears to be correct, using binary variables as Garns intended, in about half a second. By contrast, with GLPK, glpsol --mps sudoku.mps --output sudoku.sol finds a solution in 4 seconds. The CBC user’s guide explains in detail how CBC supports integer variables but does not explain how to use the cbc executable; Prof. Haroldo Gambini Santos wrote a short separate CBC command-line guide in 2011; it recommends the simpler

$ cbc sudoku.mps solve solu sudoku.solution

which seems to be equivalent. Without the “solve” command (or the lower-level dualsimplex and branch commands in the version above), cbc will still be happy to write sudoku.solution, but it may contain no solution, or only a solution of the linear relaxation of the problem.

On a more complex problem, unfortunately, I got CBC to dump core with this command. The other command worked better but I think it may not have actually given the right answer.

A nice feature of CBC is that if you kill it with ^C it responds by displaying information about the process it’s made so far on the problem (and then, as you would expect, exiting).

CBC also accepts input in CPLEX LP format.

It seems that CLP, CBC, and COIN-OR in general have stagnated somewhat since John J. Forrest, the primary author of CLP and CBC, has retired from IBM Research around 2005; there’s mojibake in the CLP User Guide and the CBC User Guide, for example. CLP and CBC are primarily intended to be invoked via an API, but they support MPS input.

CLP, and perhaps COIN-OR as a whole, descends from a proprietary IBM product called OSL, “Optimization Solutions and Library”, which had its last release in 2000 and was withdrawn in 2003 with a recommendation that users switch to COIN-OR:

For those customers considering options in the area of Optimization, IBM is strongly recommending the use of the open-source COIN solution. Much of the IBM Research Division’s efforts in Optimization for the past few years (along with effort from many others in the open-source community) have gone into developing and enhancing COIN. It has progressed to the point where IBM Research feels that COIN is equal or superior to OSL in most respects....and COIN is very competitive functionally and from a performance viewpoint with other commercially available Optimization offerings.

According to the user guide, running make unitTest builds an executable called clp which can be used as a standalone solver with MPS input. Presumably the cbc executable is similar.

(I suspect IBM acquired ILOG and thus CPLEX around 2003; the announcement above suggests that it was maybe a bit after 2003, since it doesn’t mention CPLEX.)

DSDP (solver)

The abstract of the DSDP tech report (ANL/MCS-TM-277) says:

DSDP implements the dual-scaling algorithm for semidefinite programming. The source code if[sic] this interior-point solver, written entirely in ANSI C, is freely available. The solver can be used as a subroutine library, as a function within the MATLAB environment, or as an executable that reads and writes to files. Initiated in 1997, DSDP has developed into an efficient and robust general purpose solver for semidefinite programming. Although the solver is written with semidefinite programming in mind, it can also be used for linear programming and other constraint cones.

Apparently semidefinite programming is some kind of generalization of linear programming:

All linear programs can be expressed as SDPs, and via hierarchies of SDPs the solutions of polynomial optimization problems can be approximated. ... A linear programming problem is one in which we wish to maximize or minimize a linear objective function of real variables over a polytope. In semidefinite programming, we instead use real-valued vectors and are allowed to take the dot product of vectors; nonnegativity constraints on real variables in LP (linear programming) are replaced by semidefiniteness constraints on matrix variables in SDP (semidefinite programming).

The paper that introduced semidefinite programming gives further context.

Semidefinite programming unifies several standard problems (eg, linear and quadratic programming) and finds many applications in engineering. Although semidefinite programs are much more general than linear programs, they are just as easy to solve. Most interior-point methods for linear programming have been generalized to semidefinite programs.

This all sounds great, but, to me at least, it is not at all apparent how to use the dsdp5 command to solve a linear program. The documentation is all for the library’s C API, and so are the examples.

NLopt (solver)

NLopt is a nonlinear optimization library, so it probably doesn’t really belong in this list, but in theory it should be able to solve linear optimization problems too, as a special case; it will be interesting to compare it. It has bindings in C, C++, Fortran, Octave, Python, Guile, and R.

Octave (programming environment)

Although Octave is mostly a matrix computation system, it includes support for linear programming, as well as quadratic programming, general nonlinear programming, and linear least-squares solution of matrices.

Sagemath (programming environment with modeling embedded DSL)

The Sage math programming environment has a mixed integer linear programming module which provides an embedded DSL in the Sage environment for constructing MILP optimization problems. It looks a little clumsier than MathProg. It can use GLPK, CBC, CVXOPT, PPL, or some proprietary solvers.

It also has a semidefinite programming module, using CVXOPT as the default (and by default only supported) solver.

I haven’t tried it.

Parma Polyhedra Library “PPL” (solver)

The GPL C++ Parma Polyhedra Library provides a kind of generalization of mixed integer linear programming that I don’t fully understand, especially suited for static analysis of hardware systems, with APIs in C, Java, OCaml, and Prolog. I haven’t tried it. The comprehensive user guide says a lot of things like this:

Note: The affine dimension k <= n of an NNC polyhedron P in Pn must not be confused with the space dimension n of P, which is the dimension of the enclosing vector space R^n.

CVXOPT (modeling embedded DSL and solver)

CVXOPT is a GPL3+ library for Python under active development since 2004; it is intended for convex optimization but provides a lot of general numerical functionality, including interfaces to GLPK, plus its own solvers for linear and quadratic optimization (under the rubric “cone programming”), as well as nonlinear optimization and semidefinite programming. I haven’t tried it but I assume it is suitable for writing your models in too.

Pyomo (modeling embedded DSL)

Pyomo, formerly “Coopr”, is an algebraic modeling language like ZIMPL or MathProg, but realized as an embedded DSL in Python, or perhaps more accurately, a Python API, with the concepts closely matching those of ZIMPL or MathProg. It supports linear optimization and also quadratic optimization, general nonlinear optimization, etc. It’s able to read model parameters from data files in MathProg format, but not to read MathProg models.

Pyomo can use GLPK, CBC, or some undocumented set of other solvers.

FLOPC++ (modeling embedded DSL)

FLOPC++ — another part of the COIN-OR project — is, similarly, an algebraic modeling language realized as an embedded DSL in C++. I haven’t tried it but it looks like roughly the same level of pain as Pyomo. It can at least use CBC.

Google OR-Tools (solvers and embedded DSLs in multiple languages)

Google has a set of free-software “operations research tools” called OR-Tools, including a linear optimizer called Glop, licensed under the Apache license 2.0, with bindings in Python, C++, Java, and C#. Glop doesn’t seem to support MIP, but the bundle also includes a constraint solver and a SAT solver, as well as interfaces to external MIP solvers (they recommend SCIP). I haven’t tried any of them yet.

lp_solve (solver)

lp_solve is a simplex-method solver, using branch-and-bound for mixed integer programming, which supports the MPS file format. Early versions were proprietary, but recent versions are free software.

I had some trouble getting it to interoperate with GLPK or ZIMPL. I did eventually get it to read an MPS file generated by ZIMPL, and I think I know why the GLPK-generated version was failing; my example problem (given in MathProg and LP format above) is a maximization problem, not a minimization problem. Upon translating it to MPS with zimpl -t mps min.zpl (see below), ZIMPL issued the following warning:

--- Warning: Objective function inverted to make
             minimization problem for MPS output

And then lp_solve was able to solve it with lp_solve -mps min.mps, and indeed the value of the objective function was -3.66...7 rather than 3.66...7. So I guess that the problem was that MPS doesn’t have a way to say you’re trying to maximize the objective function, not minimize it.

With this in mind, I was then able to solve the free-form MPS file I had generated with glpsol -m min.mod --wfreemps min.fmps by using lp_solve -max -fmps min.fmps, -max being an lp_solve flag to override this. GLPK can also generate fixed MPS files lp_solve can read with its -mps flag.

lp_solve has its own “LP format” which is not the same as the CPLEX LP format used by GLPK, ZIMPL, CLP, and CBC, leading to much confusion on my part. The documentation does explain this:

The lp-format is lpsolves [sic] native format to read and write lp models. Note that this format is not the same as the CPLEX LP format (see CPLEX lp files) or the Xpress LP format (see Xpress lp files)

I was able to translate the file from MPS to this incompatible LP format with lp_solve -max -fmps min.fmps -wlp min.lps.lp:

/* min */

/* Objective function */
max: +2 y +3 z;

/* Constraints */
a: +3 y +2 z <= 3;

/* Variable bounds */
y <= 4;
z <= 1;

This can be reduced to the following, a format which is quite reasonable for writing by hand, and still work as input to lp_solve:

max: +2 y +3 z;
a: +3 y +2 z <= 3;
y <= 4;
z <= 1;

Unfortunately nothing but lp_solve supports this format. (But I think it can translate it to MPS format.)

In theory lp_solve supports CPLEX LP format with the flags -rxli xli_CPLEX to read or -wxli xli_CPLEX to write, but I think the Debian package failed to include the relevant module, so this doesn’t work. (If you really needed to do this, you could use GLPK to translate from LP to MPS so lp_solve can handle it.)

lp_solve even supposedly has an “xli” to read MathProg, presumably linked with GLPK.

lp_solve seems to be much weaker than GLPK’s solvers; I have an underconstrained 9x9 Sudoku puzzle, using the above solver, here that GLPK can solve in 4.1 seconds, but lp_solve chewed on it for 22 hours without solving it before I killed it. The two projects started about the same time, but GLPK has seen continued development. I suspect that any current interest in lp_solve is largely because GLPK is GPL, so you can’t build proprietary software with it, while lp_solve is LGPL, so you can.

Like GLPK, lp_solve also has an API as well as file formats — in C, Java, and Free Pascal.

Hacks to expand what you can do with a linear optimizer

The description from lp_solve makes linear optimization sound really weak and limiting:

Solve Ax⃗ ≥ V⃗₁, with V⃗₂·x⃗ maximal.

And, in some ways, it is; but as I’m learning in this class, there are standard hacks for turning a much larger variety of problems into linear problems. Some of these are done for you automatically by model translators for algebraic modeling languages, while others aren’t.

Minimizing vs. maximizing

The simplest hack is switching between minimizing and maximizing. The lp_solve definition above might suggest that you can’t minimize things, but of course when x⃗ is chosen to make V⃗₂·x⃗ minimal, then (-1 V⃗₂) · x⃗ is maximized. Moreover, its maximal value is precisely the opposite of the minimal value for V⃗₂·x⃗; so if your solver wants to maximize, you can just flip the sign on its objective vector and the objective result, and you can minimize.

As mentioned earlier, ZIMPL will do this when generating an MPS file, since apparently MPS can’t express whether you want to maximize or minimize your objective function, and so given just an MPS file, GLPK, lp_solve, CLP, and CBC assume you want to minimize.

Mixing ≥ with ≤

Similarly, Ax⃗ ≥ V⃗₁ suggests that all your constraints must have decision variables on the left side of a “≥”. But that’s silly; the constraint 2y ≥ 4, for example, is equivalent to -2y ≤ -4. So by flipping the signs on a row of the matrix and the corresponding limit, you can effectively get both ≤ and ≥ constraints in the same matrix. And that’s how you can get criteria like 1 ≤ 2y + 3z ≤ 2.

Indeed, even the CPLEX LP format supports this directly for decision variables:

Bounds
 0 <= y <= 4
 0 <= z <= 1

Decision variables on both sides

Similarly, if you have a criterion like 2y ≤ 3z + 1, you might think that doesn’t fit into the Ax⃗ ≥ V⃗₁ mold. But all you have to do is subtract 3z from both sides — 2y - 3z ≤ 1 — and Bob’s your auntie!

Equality constraints

Similarly, if you have an equality criterion like 3p = 4q - 2r, you can convert it into two inequalities — 3p ≤ 4q - 2r and 3p ≥ 4q - 2r — which can only be simultaneously true when strict equality is satisfied. This reduces the dimensionality of your feasible region below the full dimensionality of the decision-variable space, though that’s not a problem for the simplex method.

Binary gates with big M

A common nonlinearity that can be handled fairly easily by integer linear programming is a sort of “something-or-nothing” criterion; for example, with solar panels, the amount of power you can produce at a site (and thus money you can make) is proportional to the amount of panel area installed at the site, but to produce any power at all at the site, you need to run power lines to the site, which has some fixed cost.

If you have some upper limit M on the power you can produce, you can multiply that upper limit by a binary variable W that encodes whether or not a site is active; for example, in MathProg notation:

param M {s in Sites} := available_area[s] * power_per_area;
var W {s in Sites}, binary;
subject to power_installation {s in Sites}:
    0 <= panels_installed[s] <= W[s] * M[s];

When Wₛ is 0, this forces panels_installed[s] to also be 0; when Wₛ is 1, panels_installed[s] can have any value at all, as long as it’s less than M, which was chosen to be large enough not to be a restriction in practice. (Solvers can suffer precision problems if M is too many orders of magnitude larger than the actual values of the variable, though.)

Then, you can use this variable W in your objective function to take into account the installation cost.

(I mean “gate” in the sense of something that can block something else from happening, not in the sense of a NAND gate.)

Absolute values

A much cooler trick that doesn’t even require dropping back to integer linear programming is encoding absolute values. Suppose one of the terms of an objective function you’re trying to minimize is the absolute value of some linear function f of your decision variables, |f(x⃗)|. That’s not linear, or even differentiable, because it changes direction abruptly when f(x⃗) = 0, so you can’t solve it directly with a linear solver. What you can do is replace |f(x⃗)| with a new variable t which is constrained to be greater than or equal to |f(x⃗)| using two new constraints: tf(x⃗) and t ≥ -f(x⃗). This ensures that t ≥ |f(x⃗)|, which, in conjunction with the minimization criterion, forces t = |f(x⃗)|, and again, Bob’s your auntie!

When you’re trying to maximize a sum including |f(x⃗)|, which is a more unusual thing to do, this becomes trickier; you must drop back to mixed integer programming and introduce a binary variable W and maximum M bound and use a variant of the binary-gate trick to prevent t from ever being strictly greater than f(x⃗), with t - f(x⃗) ≤ WM and t - (-f(x⃗)) ≤ (1-W)M.

Mutual exclusion

As demonstrated in the Sudoku example above, you can use a simple sum condition to require that precisely one of a set of options be chosen: W1 + W2 + W3 = 1.

Piecewise-linear functions

Mutual exclusion allows you to use piecewise-linear functions in your linear constraints and objectives. You can’t do this in the obvious way — Wf₁(x) + Wf₂(x) + Wf₃(x), with the fᵢ being the linear pieces — because it wouldn’t be linear. But you can introduce new “abscissa” variables z₁, z₂, z₃ that are gated by the Wᵢ and restricted to their own regions: 3W₁ ≤ z₁ ≤ 6W₁, for example, restricts z₁ to [3, 6] when W₁ = 1, and 0 otherwise. Then f₁(z₁) + f₂(z₂) + f₃(z₃) gives you the piecewise-linear function you wanted, and z₁ + z₂ + z₃ gives you its abscissa, which perhaps you want to constrain to be equal to x.

This is a variant of the binary-gate trick that requires no big M.

(I've written it here with single scalars zᵢ, but you can extend it to piecewise-linear functions of multiple variables in a few different ways.)

This has the disadvantage that if your piecewise-linear function is discontinuous, it ceases to be a function at the discontinuities; it can take two different values.

(I think there’s also a way to make a convex downwards, continuous piecewise-linear function out of absolute values without involving integer linear programming, in the case where you’re minimizing, but I haven’t seen it worked out yet.)

Ratios between linear forms where the denominator has known sign

Topics