How to solve LCP problem (linear complementarity problem) in python?

Is there a good library for numerically solving LCP in python?

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

+6
python
source share
4 answers

For quadratic programming with Python, I use qp -solver from cvxopt ( source ). Using this, it is easy to translate the LCP problem into the 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') 

Note:

  • the code is not fully verified, please check carefully;
  • there are certainly better solution methods than converting LCP to QP.
+4
source share

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

+3
source share

OpenOpt has a free LCP receiver written in Python + NumPy, see http://openopt.org/LCP

+1
source share

The best MCP solution algorithm (mixed nonlinear complementarity problems, more general than LCP) is the PATH solver: http://pages.cs.wisc.edu/~ferris/path.html p>

The PATH solver is available in Matlab and GAMS. Both come with a python API. I decided to use GAMS because there is a free version. So here is a step-by-step solution for a LCP solution with the python GAMS API. I used python 3.6:

  • Download and install GAMS: https://www.gams.com/download/

  • Install the API on python, as here: https://www.gams.com/latest/docs/API_PY_TUTORIAL.html I used conda, changed the directory so that there are python 3.6 applications and introduced

     python setup.py install 
  • Create a .gms file (GAMS file) lcp_for_py.gms containing:

     sets i; alias(i,j); parameters m(i,i),b(i); $gdxin lcp_input $load imb $gdxin positive variables z(i); equations Linear(i); Linear(i).. sum(j,m(i,j)*z(j)) + b(i) =g= 0; model lcp linear complementarity problem/Linear.z/; options mcp = path; solve lcp using mcp; display zL; 
  • Your python code is as follows:

     import gams #Set working directory, GamsWorkspace and the Database worDir = "<THE PATH WHERE YOU STORED YOUR .GMS-FILE>" #"C:\documents\gams\" ws=gams.GamsWorkspace(working_directory=worDir) db=ws.add_database(database_name="lcp_input") #Set the matrix and the vector of the LCP as lists matrix = [[1,1],[2,1]] vector = [0,-2] #Create the Gams Set index = [] for k in range(0,len(matrix)): index.append("i"+str(k+1)) i = db.add_set("i",1,"number of decision variables") for k in index: i.add_record(k) #Create a Gams Parameter named m and add records m = db.add_parameter_dc("m", [i,i], "matrix of the lcp") for k in range(0,len(matrix)): for l in range(0,len(matrix[0])): m.add_record([index[k],index[l]]).value = matrix[k][l] #Create a Gams Parameter named b and add records b = db.add_parameter_dc("b",[i],"bias of quadratics") for k in range(0, len(vector)): b.add_record(index[k]).value = vector[k] #run the GamsJob using the Gams File and the database lcp = ws.add_job_from_file("lcp_for_py.gms") lcp.run(databases = db) #Save the solution as a list an print it z = [] for rec in lcp.out_db["z"]: z.append(rec.level) print(z) 
+1
source share

All Articles