Come risolvere un LCP (problema di complementarietà lineare) in python?

Question

Esiste una buona libreria per risolvere numericamente un LCP in python?

Modifica: Ho bisogno di un esempio di codice Python funzionante perché la maggior parte delle librerie sembra risolvere solo problemi quadratici e ho problemi a convertire un LCP in un QP.

Answer 1

Per la programmazione quadratica con Python, utilizzo lo qp -solver da cvxopt (source). Usando questo, è semplice tradurre il problema LCP in un problema di QP (vedere Wikipedia). Esempio:

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') 

Si prega di notare:

• il codice è completamente testato, si prega di verificare con attenzione;
• ci sono sicuramente tecniche di soluzione migliori rispetto alla trasformazione di LCP in QP.

Answer 2

Dai un'occhiata allo scikit OpenOpt. Ha un esempio di fare quadratic programming e credo che vada oltre le routine di ottimizzazione di SciPy. NumPy è richiesto per utilizzare OpenOpt. Credo che la pagina di wikipedia che ci hai indicato per LCP descriva come risolvere un LCP tramite QP.

Answer 3

OpenOpt ha un LCP gratuito solver scritto in Python + NumPy vedere http://openopt.org/LCP

Answer 4

il miglior algoritmo per la soluzione MCP (problemi di complementarità non lineari misti, più generali di LCP) è il risolutore PATH: http://pages.cs.wisc.edu/~ferris/path.html

Il risolutore PATH è disponibile in MATLAB e GAMS. Entrambi vengono con un'API python. Ho scelto di utilizzare GAMS perché esiste una versione gratuita. Quindi ecco una soluzione passo passo per risolvere un LCP con l'API python di GAMS. Ho usato python 3.6:

1. scaricare e installare GAMS: https://www.gams.com/download/

2. installare l'API a Python come qui: https://www.gams.com/latest/docs/API_PY_TUTORIAL.html ho usato Conda, cambiato la directory sono stati i apifiles di Python 3.6 sono stati e sono entrato

   python setup.py install 
   

3. Creare un .gms file (file GAMS) lcp_for_py.gms contenenti:

   sets i; 
   
   alias(i,j); 
   
   parameters m(i,i),b(i); 
   
   $gdxin lcp_input 
   $load i m b 
   $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 z.L; 
   

4. Il codice Python è come questo:

   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)