Linear Programming in Python with PuLP

Update: a much better solution is to use CVXOPT. See this follow-up post for details.

In this post, we will see how to solve a Linear Program (LP) in Python. As an example, we suppose that we have a set of affine functions \(f_i({\bf x}) = a_i + {\bf b}_i^\top {\bf x}\), and we want to make all of them as small as possible, that is to say, to minimize their maximum. We will formulate this problem as an LP program, see how to solve it in Python using the PuLP library, and finally compare the performances of three LP solvers. As bonus, you will also find installation instructions for CoinMP in the end, since it is not straightforward.

PuLP setup

On Debian or Ubuntu, installing PuLP is relatively easy:

$ sudo pip install pulp            # PuLP
$ sudo apt-get install glpk-utils  # GLPK
$ sudo apt-get install coinor-cbc  # CoinOR

From there, you should be able to import pulp from Python.

Linear Problem

We can generate a random instance of our linear problem as follows:

from pylab import dot, random
n = 50
a = 2. * random(n) - 1.
B = 2. * random((n, n)) - 1.
f = lambda (i, x): a[i] + dot(B[i], x)
objective = lambda x: max([f(i, x) for i in range(n)])

The goal is now to find a vector \({\bf x}\) of length \(n\) such that \(\mathrm{objective}({\bf x})\) attains its minimum value. We turn this into an LP problem by introducing an additional variable \(m\) to represent the maximum:

\begin{equation*} \begin{array}{rl} \textrm{minimize} & m \\ \textrm{subject to} & \forall i, \ a_i + \sum_j b_{ij}\,x_j \leq m \\ & \forall i, x_\min \leq x_i \leq x_\max \\ \end{array} \end{equation*}

The objective is now to find a couple \(({\bf x}, m)\) such that all constraints are satisfied and \(m\) is as small as possible. Note that I added the constraint that \({\bf x} \in [x_\min, x_\max]^n\), just for fun. Here is a Python function that solves this LP using PuLP:

import pulp

def solve_minmax(n, a, B, x_min=None, x_max=None):
    x = pulp.LpVariable.dicts("x", range(n + 1), x_min, x_max)
    lp_prob = pulp.LpProblem("Minmax Problem", pulp.LpMinimize)
    lp_prob += pulp.lpSum([x[n]]), "Minimize_the_maximum"

    for i in range(n):
        label = "Max_constraint_%d" % i
        dot_B_x = pulp.lpSum([B[i][j] * x[j] for j in range(n)])
        condition = pulp.lpSum([x[n]]) >= a[i] + dot_B_x
        lp_prob += condition, label

    lp_prob.writeLP("MinmaxProblem.lp")  # optional
    lp_prob.solve()

    print "Status:", pulp.LpStatus[lp_prob.status]
    for v in lp_prob.variables():
        print v.name, "=", v.varValue
    print "Total Cost =", pulp.value(lp_prob.objective)

Note that we encoded the additional variable \(m\) as \({\bf x}_n\). Also, setting x_min and x_max to None will disable the lower- and upper-bound constraints.

Using various LP solvers

PuLP is actually a wrapping library that generates input files for various LP solvers, for instance GLPK or COIN. To check what solvers are available on your machine, run:

In [31]: pulp.pulpTestAll()
Solver pulp.solvers.PULP_CBC_CMD unavailable.
Solver pulp.solvers.CPLEX_DLL unavailable.
Solver pulp.solvers.CPLEX_CMD unavailable.
Solver pulp.solvers.CPLEX_PY unavailable.
    Testing zero subtraction
    Testing continuous LP solution
    Testing maximize continuous LP solution
    ...
* Solver pulp.solvers.COIN_CMD passed.
Solver pulp.solvers.COINMP_DLL unavailable.
    Testing zero subtraction
    Testing continuous LP solution
    Testing maximize continuous LP solution
    ...
* Solver pulp.solvers.GLPK_CMD passed.
Solver pulp.solvers.XPRESS unavailable.
Solver pulp.solvers.GUROBI unavailable.
Solver pulp.solvers.GUROBI_CMD unavailable.
Solver pulp.solvers.PYGLPK unavailable.
Solver pulp.solvers.YAPOSIB unavailable.

Here on my machine, only COIN_CMD and GLPK_CMD passed the tests. To use one or the other, replace the call to solve() in the function above by the solver you want to use, for instance:

lp_prob.solve(pulp.COIN_CMD())

Performance comparison

Easy access to all these solvers opens the way for testing: which one is the fastest? For this purpose, let us rewrite solve_minmax() quickly:

def solve_random_minmax(n, solver):
    x = pulp.LpVariable.dicts("x", range(n + 1), -42, +42)
    lp_prob = pulp.LpProblem("Compute internal torques", pulp.LpMinimize)
    lp_prob += pulp.lpSum([x[n]]), "Minimize_the_maximum"
    a, B = 2 * pylab.random(n) - 1., 2 * pylab.random((n, n)) - 1.
    for i in range(n):
        label = "Max_constraint_%d" % i
        dot_B_x = pulp.lpSum([B[i][j] * x[j] for j in range(n)])
        condition = pulp.lpSum([x[n]]) >= a[i] + dot_B_x
        lp_prob += condition, label
    lp_prob.solve(solver)

On my machine, I get the following performances:

In [42]: %timeit solve_random_minmax(50, pulp.COIN_CMD())
10 loops, best of 3: 113 ms per loop

In [43]: %timeit solve_random_minmax(50, pulp.GLPK_CMD())
10 loops, best of 3: 120 ms per loop

In [44]: %timeit pylab.random(50), pylab.random((50, 50))
10000 loops, best of 3: 42.6 us per loop

(The last %timeit checks that the time taken to generate random inputs is negligible.) Judging from the CMD prefix in their names, I was wondering how much time was lost communicating with the LP solvers via the command-line. I thus went through the trouble of installing CoinMP (instructions below). As it turned out, COINMP_DLL was around 10-15% faster than the CMD solvers for n=50 on my machine.

In [112]: %timeit solve_random_minmax(50, pulp.COINMP_DLL())
10 loops, best of 3: 100 ms per loop

The improvement does not seem to vanish for larger values of \(n\):

In [118]: %timeit -n10 solve_random_minmax(500, pulp.COINMP_DLL())
10 loops, best of 3: 10.9 s per loop

In [119]: %timeit -n10 solve_random_minmax(500, pulp.COIN_CMD())
10 loops, best of 3: 12 s per loop

In [120]: %timeit -n10 solve_random_minmax(500, pulp.GLPK_CMD())
10 loops, best of 3: 13.6 s per loop

Installing CoinMP and PuLP

CoinMP is not available by default when you install PuLP, and getting it to work was painful for me. Here are the installation instructions that I have figured out.

  1. Install pulp using easy_install:
$ sudo easy_install pulp
  1. Check-out the pulp-or repositoring and build it (instructions found here):
$ svn checkout http://pulp-or.googlecode.com/svn/trunk/ pulp-or
$ cd pulp-or
$ ./configure
$ make
$ make test
$ chmod +x install-sh  # wrong permissions in the repository
$ make install
  1. Copy the library files from lib/ to /usr/local/lib/python*/dist-packages/PuLP*/pulp/solverdir. On my machine, it was:
$ sudo cp ./lib/*.so
/usr/local/lib/python2.6/dist-packages/PuLP-1.5.4-py2.6.egg/pulp/solverdir/
  1. Add the following line to pulp.cfg.linux (which is in the folder above) in the "[locations]" section. Make sure you update the location for "liblapack.so" if you have installed it elsewhere on your system.
CoinMPPath = /usr/lib/liblapack.so, %(here)s/solverdir/libCoinUtils.so,
%(here)s/solverdir/libClp.so, %(here)s/solverdir/libOsi.so,
%(here)s/solverdir/libOsiClp.so, %(here)s/solverdir/libCgl.so,
%(here)s/solverdir/libCbc.so, %(here)s/solverdir/libCbcSolver.so,
%(here)s/solverdir/libCoinMP.so

And voilà :) You should be able to use the solver. However, don't expect the function pulpTestAll() to work properly for COINMP_DLL as it is not fully supported in PuLP (as of version 1.5.4). You should also be cautious that CoinMP has some known bugs that are still pending (I hear: "incorrect status messages for integer infeasible problems").

Content on this website is under the CC-BY 4.0 license.