views:

226

answers:

3

Is there a good library to numericly solve an LCP in python ?

Edit: I need a working python code example because most libraries seem to only solve quadratic problems and i have problems converting an LCP into a QP.

A: 

Take a look at numpy.

KangOl
Can you give me an example of how to convert an LCP into a QP to be able to solve it by numpy ? I need this for computing 2d contact forces in a simulation.
abenthy
numpy doesn't offer solvers
stephan
+1  A: 

Take a look at the scikit OpenOpt. It has an example of doing quadratic programming and I believe that it goes beyond SciPy's optimization routines. NumPy is required to use OpenOpt. I believe that the wikipedia page that you pointed us to for LCP describes how to solve a LCP by QP.

Justin Peel
OpenOpt is a good suggestion. Actually, it offers the cvxopt suggested by me as one of its solvers.
stephan
+2  A: 

For quadratic programming with Python, I use the qp-solver from cvxopt. Using this, it is straightforward to translate the LCP problem into a QP problem (see Wikipedia). Example:

from cvxopt import matrix, spmatrix
from cvxopt.blas import gemv
from cvxopt.solvers import qp

def append_matrix_at_bottom(A, B):
    l = []
    for x in xrange(A.size[1]):
        for i in xrange(A.size[0]):
            l.append(A[i+x*A.size[0]])
        for i in xrange(B.size[0]):
            l.append(B[i+x*B.size[0]])
    return matrix(l,(A.size[0]+B.size[0],A.size[1]))

M = matrix([[ 4.0,  6, -4,    1.0 ],
            [ 6,  1,  1.0,     2.0 ],
            [-4,  1.0,   2.5,  -2.0 ],
            [ 1.0,   2.0,   -2.0,     1.0 ]])
q = matrix([12, -10, -7.0, 3])

I = spmatrix(1.0, range(M.size[0]), range(M.size[1]))
G = append_matrix_at_bottom(-M, -I)   # inequality constraint G z <= h
h = matrix([x for x in q] + [0.0 for _x in range(M.size[0])])

sol = qp(2.0 * M, q, G, h)      # find z, w, so that w = M z + q
if sol['status'] == 'optimal':
    z = sol['x']
    w = matrix(q)
    gemv(M, z, w, alpha=1.0, beta=1.0)   # w = M z + q
    print(z)
    print(w)
else:
    print('failed')

Please note:

  • the code is totally untested, please check carefully;
  • there surely are better solution techniques than transforming LCP into QP.
stephan