Detailed Description

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] (see pp. 224-230 in the latter) 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.

#include <drake/solvers/moby_lcp_solver.h>

Public Member Functions
	MobyLcpSolver ()
	~MobyLcpSolver () final
bool	SolveLcpFast (const Eigen::MatrixXd &M, const Eigen::VectorXd &q, Eigen::VectorXd *z, double zero_tol=-1) const
	Fast pivoting algorithm for LCPs of the form M = PAPᵀ, q = Pb, where b ∈ ℝᵐ, P ∈ ℝⁿˣᵐ, and A ∈ ℝᵐˣᵐ (where A is positive definite).
bool	SolveLcpFastRegularized (const Eigen::MatrixXd &M, const Eigen::VectorXd &q, Eigen::VectorXd *z, int min_exp=-20, unsigned step_exp=4, int max_exp=20, double zero_tol=-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).
bool	SolveLcpLemke (const Eigen::MatrixXd &M, const Eigen::VectorXd &q, Eigen::VectorXd *z, double piv_tol=-1, double zero_tol=-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.
bool	SolveLcpLemkeRegularized (const Eigen::MatrixXd &M, const Eigen::VectorXd &q, Eigen::VectorXd *z, int min_exp=-20, unsigned step_exp=1, int max_exp=1, double piv_tol=-1, double zero_tol=-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.
int	get_num_pivots () const
	Returns the number of pivoting operations made by the last LCP solve.
void	reset_num_pivots ()
	Resets the number of pivoting operations made by the last LCP solver to zero.
Does not allow copy, move, or assignment
	MobyLcpSolver (const MobyLcpSolver &)=delete
MobyLcpSolver &	operator= (const MobyLcpSolver &)=delete
	MobyLcpSolver (MobyLcpSolver &&)=delete
MobyLcpSolver &	operator= (MobyLcpSolver &&)=delete
Public Member Functions inherited from SolverBase
	~SolverBase () override
MathematicalProgramResult	Solve (const MathematicalProgram &prog, const std::optional< Eigen::VectorXd > &initial_guess=std::nullopt, const std::optional< SolverOptions > &solver_options=std::nullopt) const
	Like SolverInterface::Solve(), but the result is a return value instead of an output argument.
void	Solve (const MathematicalProgram &, const std::optional< Eigen::VectorXd > &, const std::optional< SolverOptions > &, MathematicalProgramResult *) const override
	Solves an optimization program with optional initial guess and solver options.
bool	available () const override
	Returns true iff support for this solver has been compiled into Drake.
bool	enabled () const override
	Returns true iff this solver is properly configured for use at runtime.
SolverId	solver_id () const final
	Returns the identifier of this solver.
bool	AreProgramAttributesSatisfied (const MathematicalProgram &) const override
	Returns true iff the program's attributes are compatible with this solver's capabilities.
std::string	ExplainUnsatisfiedProgramAttributes (const MathematicalProgram &) const override
	Describes the reasons (if any) why the program is incompatible with this solver's capabilities.
	SolverBase (const SolverBase &)=delete
SolverBase &	operator= (const SolverBase &)=delete
	SolverBase (SolverBase &&)=delete
SolverBase &	operator= (SolverBase &&)=delete
Public Member Functions inherited from SolverInterface
virtual	~SolverInterface ()
	SolverInterface (const SolverInterface &)=delete
SolverInterface &	operator= (const SolverInterface &)=delete
	SolverInterface (SolverInterface &&)=delete
SolverInterface &	operator= (SolverInterface &&)=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.

Static versions of the instance methods with similar names.
MathematicalProgramResult	Solve (const MathematicalProgram &prog, const std::optional< Eigen::VectorXd > &initial_guess=std::nullopt, const std::optional< SolverOptions > &solver_options=std::nullopt) const
	Like SolverInterface::Solve(), but the result is a return value instead of an output argument.
void	Solve (const MathematicalProgram &, const std::optional< Eigen::VectorXd > &, const std::optional< SolverOptions > &, MathematicalProgramResult *) const override
	Solves an optimization program with optional initial guess and solver options.
static SolverId	id ()
static bool	is_available ()
static bool	is_enabled ()
static bool	ProgramAttributesSatisfied (const MathematicalProgram &)

Additional Inherited Members
Protected Member Functions inherited from SolverBase
	SolverBase (const SolverId &id, std::function< bool()> available, std::function< bool()> enabled, std::function< bool(const MathematicalProgram &)> are_satisfied, std::function< std::string(const MathematicalProgram &)> explain_unsatisfied=nullptr)
	Constructs a SolverBase with the given default implementations of the solver_id(), available(), enabled(), AreProgramAttributesSatisfied(), and ExplainUnsatisfiedProgramAttributes() methods.
virtual void	DoSolve2 (const MathematicalProgram &prog, const Eigen::VectorXd &initial_guess, internal::SpecificOptions options, MathematicalProgramResult result) const
	(Internal use only) Like DoSolve() but using SpecificOptions instead of SolverOptions.
Protected Member Functions inherited from SolverInterface
	SolverInterface ()

Constructor & Destructor Documentation

◆ MobyLcpSolver() [1/3]

MobyLcpSolver ( const MobyLcpSolver & )

delete

◆ MobyLcpSolver() [2/3]

MobyLcpSolver ( MobyLcpSolver && )

delete

◆ MobyLcpSolver() [3/3]

MobyLcpSolver ( )

◆ ~MobyLcpSolver()

~MobyLcpSolver ( )

final

Member Function Documentation

◆ ComputeZeroTolerance()

template<class U>

U ComputeZeroTolerance ( const MatrixX< U > & M )

static

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

◆ get_num_pivots()

int get_num_pivots ( ) const

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

◆ id()

SolverId id ( )

static

◆ is_available()

bool is_available ( )

static

◆ is_enabled()

bool is_enabled ( )

static

◆ operator=() [1/2]

MobyLcpSolver & operator= ( const MobyLcpSolver & )

delete

◆ operator=() [2/2]

MobyLcpSolver & operator= ( MobyLcpSolver && )

delete

◆ ProgramAttributesSatisfied()

bool ProgramAttributesSatisfied ( const MathematicalProgram & )

static

◆ reset_num_pivots()

void reset_num_pivots ( )

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

◆ Solve() [1/2]

void Solve	(	const MathematicalProgram &	prog,
		const std::optional< Eigen::VectorXd > &	initial_guess,
		const std::optional< SolverOptions > &	solver_options,
		MathematicalProgramResult *	result ) const

overridevirtual

Solves an optimization program with optional initial guess and solver options.

Note that these initial guess and solver options are not written to prog. If the prog has set an initial guess, and initial_guess is set, then initial_guess takes priority. If the prog has set an option for a solver, and solver_options contains a different value for the same option on the same solver, then solver_options takes priority. Derived implementations of this interface may elect to throw std::exception for badly formed programs.

Implements SolverInterface.

◆ Solve() [2/2]

MathematicalProgramResult Solve	(	const MathematicalProgram &	prog,
		const std::optional< Eigen::VectorXd > &	initial_guess = std::nullopt,
		const std::optional< SolverOptions > &	solver_options = std::nullopt ) const

Like SolverInterface::Solve(), but the result is a return value instead of an output argument.

◆ SolveLcpFast()

bool SolveLcpFast	(	const Eigen::MatrixXd &	M,
		const Eigen::VectorXd &	q,
		Eigen::VectorXd *	z,
		double	zero_tol = -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]	M	the LCP matrix.
[in]	q	the LCP vector.
[in,out]	z	the 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_tol	The tolerance for testing against zero. If the tolerance is negative (default) the solver will determine a generally reasonable tolerance.

Exceptions

std::exception if 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.

◆ SolveLcpFastRegularized()

bool SolveLcpFastRegularized	(	const Eigen::MatrixXd &	M,
		const Eigen::VectorXd &	q,
		Eigen::VectorXd *	z,
		int	min_exp = -20,
		unsigned	step_exp = 4,
		int	max_exp = 20,
		double	zero_tol = -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]	M	the LCP matrix.
[in]	q	the LCP vector.
[in,out]	z	the 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_exp	The 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_exp	The 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_exp	The 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_tol	The tolerance for testing against zero. If the tolerance is negative (default) the solver will determine a generally reasonable tolerance.

Exceptions

std::exception if 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.

◆ SolveLcpLemke()

bool SolveLcpLemke	(	const Eigen::MatrixXd &	M,
		const Eigen::VectorXd &	q,
		Eigen::VectorXd *	z,
		double	piv_tol = -1,
		double	zero_tol = -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]	M	the LCP matrix.
[in]	q	the LCP vector.
[in,out]	z	the 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_tol	The tolerance for testing against zero. If the tolerance is negative (default) the solver will determine a generally reasonable tolerance.
[in]	piv_tol	The 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::exception if M is not square or the dimensions of M do not match the length of q.

[Cottle 1992] R. Cottle, J.-S. Pang, and R. Stone. The Linear Complementarity Problem. Academic Press, 1992.
[LEMKE] P. Fackler and M. Miranda. LEMKE. http://people.sc.fsu.edu/~burkardt/m_src/lemke/lemke.m

◆ SolveLcpLemkeRegularized()

bool SolveLcpLemkeRegularized	(	const Eigen::MatrixXd &	M,
		const Eigen::VectorXd &	q,
		Eigen::VectorXd *	z,
		int	min_exp = -20,
		unsigned	step_exp = 1,
		int	max_exp = 1,
		double	piv_tol = -1,
		double	zero_tol = -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] z the 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.

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

drake/solvers/moby_lcp_solver.h

Detailed Description

Public Member Functions

Static Public Member Functions

Static versions of the instance methods with similar names.

Additional Inherited Members

Constructor & Destructor Documentation

◆ MobyLcpSolver() [1/3]

◆ MobyLcpSolver() [2/3]

◆ MobyLcpSolver() [3/3]

◆ ~MobyLcpSolver()

Member Function Documentation

◆ ComputeZeroTolerance()

◆ get_num_pivots()

◆ id()

◆ is_available()

◆ is_enabled()

◆ operator=() [1/2]

◆ operator=() [2/2]

◆ ProgramAttributesSatisfied()

◆ reset_num_pivots()

◆ Solve() [1/2]

◆ Solve() [2/2]

◆ SolveLcpFast()

◆ SolveLcpFastRegularized()

◆ SolveLcpLemke()

◆ SolveLcpLemkeRegularized()