# 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:

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.

- Install pulp using easy_install:

$ sudo easy_install pulp

- 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

- 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/

- 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").

## To go further

You can try out the lpsolvers module to solve linear programs with CVXOPT or other solvers available in Python (for now PuLP is not included).