Drake
MobyLCPSolver< T > Class Template Reference

A class for solving Linear Complementarity Problems (LCPs). More...

#include <drake/solvers/moby_lcp_solver.h>

Inheritance diagram for MobyLCPSolver< T >:
[legend]
Collaboration diagram for MobyLCPSolver< T >:
[legend]

Public Member Functions

 MobyLCPSolver ()=default
 
 ~MobyLCPSolver () override=default
 
void SetLoggingEnabled (bool enabled)
 
bool SolveLcpFast (const MatrixX< T > &M, const VectorX< T > &q, VectorX< T > *z, const T &zero_tol=T(-1)) const
 Fast pivoting algorithm for LCPs of the form M = PAPᵀ, q = Pb, where b ∈ ℝᵐ, P ∈ ℝⁿˣᵐ, and A ∈ ℝᵐˣᵐ (where A is positive definite). More...
 
bool SolveLcpFastRegularized (const MatrixX< T > &M, const VectorX< T > &q, VectorX< T > *z, int min_exp=-20, unsigned step_exp=4, int max_exp=20, const T &zero_tol=T(-1)) const
 Regularized version of the fast pivoting algorithm for LCPs of the form M = PAPᵀ, q = Pb, where b ∈ ℝᵐ, P ∈ ℝⁿˣᵐ, and A ∈ ℝᵐˣᵐ (where A is positive definite). More...
 
bool SolveLcpLemke (const MatrixX< T > &M, const VectorX< T > &q, VectorX< T > *z, const T &piv_tol=T(-1), const T &zero_tol=T(-1)) const
 Lemke's Algorithm for solving LCPs in the matrix class E, which contains all strictly semimonotone matrices, all P-matrices, and all strictly copositive matrices. More...
 
bool SolveLcpLemkeRegularized (const MatrixX< T > &M, const VectorX< T > &q, VectorX< T > *z, int min_exp=-20, unsigned step_exp=1, int max_exp=1, const T &piv_tol=T(-1), const T &zero_tol=T(-1)) const
 Lemke's Algorithm for solving LCPs in the matrix class E, which contains all strictly semimonotone matrices, all P-matrices, and all strictly copositive matrices. More...
 
bool SolveLcpLemke (const Eigen::SparseMatrix< double > &M, const Eigen::VectorXd &q, Eigen::VectorXd *z, double piv_tol=-1.0, double zero_tol=-1.0) const
 Lemke's Algorithm for solving LCPs in the matrix class E, which contains all strictly semimonotone matrices, all P-matrices, and all strictly copositive matrices, for the special case of sparse matrices. More...
 
bool SolveLcpLemkeRegularized (const Eigen::SparseMatrix< double > &M, const Eigen::VectorXd &q, Eigen::VectorXd *z, int min_exp=-20, unsigned step_exp=4, int max_exp=20, double piv_tol=-1.0, double zero_tol=-1.0) const
 Regularized wrapper around Lemke's Algorithm for solving LCPs in the matrix class E. More...
 
bool available () const override
 Returns true iff this solver was enabled at compile-time. More...
 
SolutionResult Solve (MathematicalProgram &prog) const override
 Sets values for the decision variables on the given MathematicalProgram prog, or: More...
 
SolverId solver_id () const override
 Returns the identifier of this solver. More...
 
int get_num_pivots () const
 Returns the number of pivoting operations made by the last LCP solve. More...
 
void reset_num_pivots ()
 Resets the number of pivoting operations made by the last LCP solver to zero. More...
 
template<>
SolutionResult Solve (MathematicalProgram &) const
 Sets values for the decision variables on the given MathematicalProgram prog, or: More...
 
Does not allow copy, move, or assignment
 MobyLCPSolver (const MobyLCPSolver &)=delete
 
MobyLCPSolveroperator= (const MobyLCPSolver &)=delete
 
 MobyLCPSolver (MobyLCPSolver &&)=delete
 
MobyLCPSolveroperator= (MobyLCPSolver &&)=delete
 
- Public Member Functions inherited from MathematicalProgramSolverInterface
 MathematicalProgramSolverInterface ()=default
 
virtual ~MathematicalProgramSolverInterface ()=default
 
 MathematicalProgramSolverInterface (const MathematicalProgramSolverInterface &)=delete
 
MathematicalProgramSolverInterfaceoperator= (const MathematicalProgramSolverInterface &)=delete
 
 MathematicalProgramSolverInterface (MathematicalProgramSolverInterface &&)=delete
 
MathematicalProgramSolverInterfaceoperator= (MathematicalProgramSolverInterface &&)=delete
 

Static Public Member Functions

template<class U >
static U ComputeZeroTolerance (const MatrixX< U > &M)
 Calculates the zero tolerance that the solver would compute if the user does not specify a tolerance. More...
 

Detailed Description

template<class T>
class drake::solvers::MobyLCPSolver< T >

A class for solving Linear Complementarity Problems (LCPs).

Solving a LCP requires finding a solution to the problem:

Mz + q = w
z ≥ 0
w ≥ 0
zᵀw = 0

(where M ∈ ℝⁿˣⁿ and q ∈ ℝⁿ are problem inputs and z ∈ ℝⁿ and w ∈ ℝⁿ are unknown vectors) or correctly reporting that such a solution does not exist. In spite of their linear structure, solving LCPs is NP-Hard [Cottle 1992]. However, some LCPs are significantly easier to solve. For instance, it can be seen that the LCP is solvable in worst-case polynomial time for the case of symmetric positive-semi-definite M by formulating it as the following convex quadratic program:

minimize:   f(z) = zᵀw = zᵀ(Mz + q)
subject to: z ≥ 0
            Mz + q ≥ 0

Note that this quadratic program's (QP) objective function at the minimum z cannot be less than zero, and the LCP is only solved if the objective function at the minimum is equal to zero. Since the seminal result of Karmarkar, it has been known that convex QPs are solvable in polynomial time [Karmarkar 1984].

The difficulty of solving an LCP is characterized by the properties of its particular matrix, namely the class of matrices it belongs to. Classes include, for example, Q, Q₀, P, P₀, copositive, and Z matrices. [Cottle 1992] and Murty 1998 describe relevant matrix classes in more detail.

  • [Cottle 1992] R. Cottle, J.-S. Pang, and R. Stone. The Linear Complementarity Problem. Academic Press, 1992.
  • [Karmarkar 1984] N. Karmarkar. A New Polynomial-Time Algorithm for Linear Programming. Combinatorica, 4(4), pp. 373-395.
  • [Murty 1988] K. Murty. Linear Complementarity, Linear and Nonlinear Programming. Heldermann Verlag, 1988.

Constructor & Destructor Documentation

MobyLCPSolver ( const MobyLCPSolver< T > &  )
delete
MobyLCPSolver ( MobyLCPSolver< T > &&  )
delete
MobyLCPSolver ( )
default
~MobyLCPSolver ( )
overridedefault

Member Function Documentation

bool available ( ) const
inlineoverridevirtual

Returns true iff this solver was enabled at compile-time.

Implements MathematicalProgramSolverInterface.

static U ComputeZeroTolerance ( const MatrixX< U > &  M)
inlinestatic

Calculates the zero tolerance that the solver would compute if the user does not specify a tolerance.

int get_num_pivots ( ) const
inline

Returns the number of pivoting operations made by the last LCP solve.

MobyLCPSolver& operator= ( MobyLCPSolver< T > &&  )
delete
MobyLCPSolver& operator= ( const MobyLCPSolver< T > &  )
delete
void reset_num_pivots ( )
inline

Resets the number of pivoting operations made by the last LCP solver to zero.

void SetLoggingEnabled ( bool  enabled)
SolutionResult Solve ( MathematicalProgram prog) const
virtual

Sets values for the decision variables on the given MathematicalProgram prog, or:

  • If no solver is available, throws std::runtime_error
  • If the solver returns an error, returns a nonzero SolutionResult.

Implements MathematicalProgramSolverInterface.

SolutionResult Solve ( MathematicalProgram prog) const
overridevirtual

Sets values for the decision variables on the given MathematicalProgram prog, or:

  • If no solver is available, throws std::runtime_error
  • If the solver returns an error, returns a nonzero SolutionResult.

Implements MathematicalProgramSolverInterface.

Here is the call graph for this function:

bool SolveLcpFast ( const MatrixX< T > &  M,
const VectorX< T > &  q,
VectorX< T > *  z,
const T &  zero_tol = T(-1) 
) const

Fast pivoting algorithm for LCPs of the form M = PAPᵀ, q = Pb, where b ∈ ℝᵐ, P ∈ ℝⁿˣᵐ, and A ∈ ℝᵐˣᵐ (where A is positive definite).

Therefore, q is in the range of P and M is positive semi-definite. An LCP of this form is also guaranteed to have a solution [Cottle 1992].

This particular implementation focuses on the case where the solution requires few nonzero nonbasic variables, meaning that few z variables need be nonzero to find a solution to Mz + q = w. This algorithm, which is based off of Dantzig's Principle Pivoting Method I [Cottle 1992] is described in [Drumwright 2015]. This algorithm is able to use "warm" starting- a solution to a "nearby" LCP can be used to find the solution to a given LCP more quickly.

Although this solver is theoretically guaranteed to give a solution to the LCPs described above, accumulated floating point error from pivoting operations could cause the solver to fail. Additionally, the solver can be applied with some success to problems outside of its guaranteed matrix class. For these reasons, the solver returns a flag indicating success/failure.

Parameters
[in]Mthe LCP matrix.
[in]qthe LCP vector.
[in,out]zthe solution to the LCP on return (if the solver succeeds). If the length of z is equal to the length of q, the solver will attempt to use z's value as a starting solution. If the solver fails (returns false), z will be set to the zero vector.
[in]zero_tolThe tolerance for testing against zero. If the tolerance is negative (default) the solver will determine a generally reasonable tolerance.
Exceptions
std::logic_errorif M is non-square or M's dimensions do not equal q's dimension.
Returns
true if the solver succeeded and false otherwise.
  • [Cottle 1992] R. Cottle, J.-S. Pang, and R. Stone. The Linear Complementarity Problem. Academic Press, 1992.
  • [Drumwright 2015] E. Drumwright. Rapidly computable viscous friction and no-slip rigid contact models. arXiv: 1504.00719v1. 2015.

Here is the call graph for this function:

bool SolveLcpFastRegularized ( const MatrixX< T > &  M,
const VectorX< T > &  q,
VectorX< T > *  z,
int  min_exp = -20,
unsigned  step_exp = 4,
int  max_exp = 20,
const T &  zero_tol = T(-1) 
) const

Regularized version of the fast pivoting algorithm for LCPs of the form M = PAPᵀ, q = Pb, where b ∈ ℝᵐ, P ∈ ℝⁿˣᵐ, and A ∈ ℝᵐˣᵐ (where A is positive definite).

Therefore, q is in the range of P and M is positive semi-definite. Please see SolveLcpFast() for more documentation about the particular algorithm.

This implementation wraps that algorithm with a Tikhonov-type regularization approach. Specifically, this implementation repeatedly attempts to solve the LCP:

(M + Iα)z + q = w
z ≥ 0
w ≥ 0
zᵀw = 0

where I is the identity matrix and α ≪ 1, using geometrically increasing values of α, until the LCP is solved. Cottle et al. describe how, for sufficiently large α, the LCP will always be solvable [Cottle 1992], p. 493.

Although this solver is theoretically guaranteed to give a solution to the LCPs described above, accumulated floating point error from pivoting operations could cause the solver to fail. Additionally, the solver can be applied with some success to problems outside of its guaranteed matrix class. For these reasons, the solver returns a flag indicating success/failure.

Parameters
[in]Mthe LCP matrix.
[in]qthe LCP vector.
[in,out]zthe solution to the LCP on return (if the solver succeeds). If the length of z is equal to the length of q, the solver will attempt to use z's value as a starting solution.
[in]min_expThe minimum exponent for computing α over [10ᵝ, 10ᵞ] in steps of 10ᵟ, where β is the minimum exponent, γ is the maximum exponent, and δ is the stepping exponent.
[in]step_expThe stepping exponent for computing α over [10ᵝ, 10ᵞ] in steps of 10ᵟ, where β is the minimum exponent, γ is the maximum exponent, and δ is the stepping exponent.
[in]max_expThe maximum exponent for computing α over [10ᵝ, 10ᵞ] in steps of 10ᵟ, where β is the minimum exponent, γ is the maximum exponent, and δ is the stepping exponent.
[in]zero_tolThe tolerance for testing against zero. If the tolerance is negative (default) the solver will determine a generally reasonable tolerance.
Exceptions
std::logic_errorif M is non-square or M's dimensions do not equal q's dimension.
Returns
true if the solver succeeded and false if the solver did not find a solution for α = 10ᵞ.
See also
SolveLcpFast()
  • [Cottle, 1992] R. Cottle, J.-S. Pang, and R. Stone. The Linear Complementarity Problem. Academic Press, 1992.

Here is the call graph for this function:

bool SolveLcpLemke ( const MatrixX< T > &  M,
const VectorX< T > &  q,
VectorX< T > *  z,
const T &  piv_tol = T(-1),
const T &  zero_tol = T(-1) 
) const

Lemke's Algorithm for solving LCPs in the matrix class E, which contains all strictly semimonotone matrices, all P-matrices, and all strictly copositive matrices.

Lemke's Algorithm is described in [Cottle 1992], Section 4.4. This implementation was adapted from the LEMKE Library [LEMKE] for Matlab; this particular implementation fixes a bug in LEMKE that could occur when multiple indices passed the minimum ratio test.

Although this solver is theoretically guaranteed to give a solution to the LCPs described above, accumulated floating point error from pivoting operations could cause the solver to fail. Additionally, the solver can be applied with some success to problems outside of its guaranteed matrix classes. For these reasons, the solver returns a flag indicating success/failure.

Parameters
[in]Mthe LCP matrix.
[in]qthe LCP vector.
[in,out]zthe solution to the LCP on return (if the solver succeeds). If the length of z is equal to the length of q, the solver will attempt to use z's value as a starting solution. This warmstarting is generally not recommended: it has a predisposition to lead to a failing pivoting sequence. If the solver fails (returns false), z will be set to the zero vector.
[in]zero_tolThe tolerance for testing against zero. If the tolerance is negative (default) the solver will determine a generally reasonable tolerance.
[in]piv_tolThe tolerance for testing against zero, specifically used for the purpose of finding variables for pivoting. If the tolerance is negative (default) the solver will determine a generally reasonable tolerance.
Returns
true if the solver believes it has computed a solution (which it determines by the ability to "pivot out" the "artificial" variable (see [Cottle 1992]) and false otherwise.
Warning
The caller should verify that the algorithm has solved the LCP to the desired tolerances on returns indicating success.
Exceptions
std::logic_errorif M is not square or the dimensions of M do not match the length of q.

Here is the call graph for this function:

bool SolveLcpLemke ( const Eigen::SparseMatrix< double > &  M,
const Eigen::VectorXd &  q,
Eigen::VectorXd *  z,
double  piv_tol = -1.0,
double  zero_tol = -1.0 
) const

Lemke's Algorithm for solving LCPs in the matrix class E, which contains all strictly semimonotone matrices, all P-matrices, and all strictly copositive matrices, for the special case of sparse matrices.

See the non-sparse version of SolveLcpLemke() for descriptions of the calling and return parameters.

Note
This function is not templatized because the pivoting operations make single-precision floating point solves untenable and because the underlying sparse linear system solver does not support AutoDiff.

Here is the call graph for this function:

bool SolveLcpLemkeRegularized ( const MatrixX< T > &  M,
const VectorX< T > &  q,
VectorX< T > *  z,
int  min_exp = -20,
unsigned  step_exp = 1,
int  max_exp = 1,
const T &  piv_tol = T(-1),
const T &  zero_tol = T(-1) 
) const

Lemke's Algorithm for solving LCPs in the matrix class E, which contains all strictly semimonotone matrices, all P-matrices, and all strictly copositive matrices.

Lemke's Algorithm is described in [Cottle 1992], Section 4.4.

This implementation wraps that algorithm with a Tikhonov-type regularization approach. Specifically, this implementation repeatedly attempts to solve the LCP:

(M + Iα)z + q = w
z ≥ 0
w ≥ 0
zᵀw = 0

where I is the identity matrix and α ≪ 1, using geometrically increasing values of α, until the LCP is solved. See SolveLcpFastRegularized() for description of the regularization process and the function parameters, which are identical. See SolveLcpLemke() for a description of Lemke's Algorithm. See SolveLcpFastRegularized() for a description of all calling parameters other than z, which apply equally well to this function.

Parameters
[in,out]zthe solution to the LCP on return (if the solver succeeds). If the length of z is equal to the length of q, the solver will attempt to use z's value as a starting solution. This warmstarting is generally not recommended: it has a predisposition to lead to a failing pivoting sequence.
See also
SolveLcpFastRegularized()
SolveLcpLemke()
  • [Cottle 1992] R. Cottle, J.-S. Pang, and R. Stone. The Linear Complementarity Problem. Academic Press, 1992.

Here is the call graph for this function:

bool SolveLcpLemkeRegularized ( const Eigen::SparseMatrix< double > &  M,
const Eigen::VectorXd &  q,
Eigen::VectorXd *  z,
int  min_exp = -20,
unsigned  step_exp = 4,
int  max_exp = 20,
double  piv_tol = -1.0,
double  zero_tol = -1.0 
) const

Regularized wrapper around Lemke's Algorithm for solving LCPs in the matrix class E.

See the non-sparse version of SolveLcpLemkeRegularized() for descriptions of the calling and return parameters.

Note
This function is not templatized because the pivoting operations make single-precision floating point solves untenable and because the underlying sparse linear system solver does not support AutoDiff.

Here is the call graph for this function:

SolverId solver_id ( ) const
overridevirtual

Returns the identifier of this solver.

Implements MathematicalProgramSolverInterface.

Here is the call graph for this function:


The documentation for this class was generated from the following files: