I have been using lpsolve in Python for a very long time, with generally good results. I have even written my own Cython wrapper to it to overcome the mess that the original Python wrapper is.
My Cython wrapper is much faster in the setup of the problem, but of course the solution time only depends on the lpsolve C code and has got nothing to do with Python.
I am only solving real-valued linear problems, no MIPs. The latest (and one of the biggest) LP I had to solve had a constraint matrix of size about 5,000 x 3,000, and the setup plus solve takes about 150 ms. The problem is, I have to solve a slightly modified version of the same problem many times in my simulation (constraints, RHS, bounds and so on are time-dependent for a simulation with many timesteps). The constraint matrix is usually very sparse, with about 0.1% - 0.5% of NNZ or less.
Using lpsolve, the setup and solve of the problem is as easy as the following:
import numpy
from lp_solve.lp_solve import PRESOLVE_COLS, PRESOLVE_ROWS, PRESOLVE_LINDEP, PRESOLVE_NONE, SCALE_DYNUPDATE, lpsolve
# The constraint_matrix id a 2D NumPy array and it contains
# both equality and inequality constraints
# And I need it to be single precision floating point, like all
# other NumPy arrays from here onwards
m, n = constraint_matrix.shape
n_le = len(inequality_constraints)
n_e = len(equality_constraints)
# Setup RHS vector
b = numpy.zeros((m, ), dtype=numpy.float32)
# Assign equality and inequality constraints (= and <=)
b[0:n_e] = equality_constraints
b[-n_le:] = inequality_constraints
# Tell lpsolve which rows are equalities and which are inequalities
e = numpy.asarray(['LE']*m)
e[0:-n_le] = 'EQ'
# Make the LP
lp = lpsolve('make_lp', m, n)
# Some options for scaling the problem
lpsolve('set_scaling', lp, SCALE_DYNUPDATE)
lpsolve('set_verbose', lp, 'IMPORTANT')
# Use presolve as it is much faster
lpsolve('set_presolve', lp, PRESOLVE_COLS | PRESOLVE_ROWS | PRESOLVE_LINDEP)
# I only care about maximization
lpsolve('set_sense', lp, True)
# Set the objective function of the problem
lpsolve('set_obj_fn', lp, objective_function)
lpsolve('set_mat', lp, constraint_matrix)
# Tell lpsolve about the RHS
lpsolve('set_rh_vec', lp, b)
# Set the constraint type (equality or inequality)
lpsolve('set_constr_type', lp, e)
# Set upper bounds for variables - lower bounds are automatically 0
lpsolve('set_upbo', lp, ub_values)
# Solve the problem
out = lpsolve('solve', lp)
# Retrieve the solution for all the variables
vars_sol = numpy.asarray([lpsolve('get_var_primalresult', lp, i) for i in xrange(m + 1, m + n + 1)], dtype=numpy.float32)
# Delete the problem, timestep done
lpsolve('delete_lp', lp)
For reasons that are too long to explain, my NumPy arrays are all single precision floating point arrays, and I'd like them to stay that way.
Now, after all this painful introduction, I would like to ask: does anyone know of another library (with reasonable Python wrappers) that allows me to setup and solve a problem of this size as fast (or potentially faster) than lpsolve? Most of the libraries I have looked at (PuLP, CyLP, PyGLPK and so on) do not seem to have a straightforward way to say "this is my entire constraint matrix, set it in one go". They mostly seem to be oriented towards being "modelling languages" that allow fancy things like this (CyLP example):
# Add variables
x = s.addVariable('x', 3)
y = s.addVariable('y', 2)
# Create coefficients and bounds
A = np.matrix([[1., 2., 0],[1., 0, 1.]])
B = np.matrix([[1., 0, 0], [0, 0, 1.]])
D = np.matrix([[1., 2.],[0, 1]])
a = CyLPArray([5, 2.5])
b = CyLPArray([4.2, 3])
x_u= CyLPArray([2., 3.5])
# Add constraints
s += A * x <= a
s += 2 <= B * x + D * y <= b
s += y >= 0
s += 1.1 <= x[1:3] <= x_u
I honestly don't care about the flexibility, I just need raw speed in the problem setup and solve. Creating NumPy matrices plus doing all those fancy operations above is definitely going to be a performance killer.
I'd rather stay on Open Source solvers if possible, but any suggestion is most welcome. My apologies for the long post.
Andrea.
