The solver package

The solver package#

QP and pQP#

class mpt4py.solvers.qp.PqpSolution(status: str, regions: mpt4py.geometry.union.polyunion.PolyUnion)[source]#

Bases: object

regions: PolyUnion#
status: str#
mpt4py.solvers.qp.QP#

alias of QuadraticProgram

class mpt4py.solvers.qp.QuadraticProgram(H: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, f: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, c: float = <factory>, pY: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, pF: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, pC: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, A: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, b: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, pB: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, Ae: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, be: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, pE: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, lb: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, ub: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, Ath: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, bth: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>)[source]#

Bases: object

Defines a multi-parametric Quadratic Program (QP) of the form:

\[\begin{split}\begin{aligned} \min\limits_x &\quad \frac{1}{2} x^\top Hx + (F \theta+f)^\top x + \theta^\top Y \theta + C\theta + c, \\ \text{s.t.} &\quad A x \leq b + B \theta, \\ & \quad A_e x = b_e + E \theta, \\ & \quad x_{lb} \leq x \leq x_{ub}, \end{aligned}\end{split}\]

where \(H\) and \(Y\) are positive semi-definite.

The upper and lower bounds are not supported in the parametric case yet.

A: ndarray[tuple[Any, ...], dtype[float64]]#
Ae: ndarray[tuple[Any, ...], dtype[float64]]#
Ath: ndarray[tuple[Any, ...], dtype[float64]]#
H: ndarray[tuple[Any, ...], dtype[float64]]#
b: ndarray[tuple[Any, ...], dtype[float64]]#
be: ndarray[tuple[Any, ...], dtype[float64]]#
bth: ndarray[tuple[Any, ...], dtype[float64]]#
c: float#
d: int#
f: ndarray[tuple[Any, ...], dtype[float64]]#
property is_parametric: bool#

Determine if the quadratic program is parametric or not.

lb: ndarray[tuple[Any, ...], dtype[float64]]#
m: int#
me: int#
n: int#
property num_eq: int#
property num_ineq: int#
pB: ndarray[tuple[Any, ...], dtype[float64]]#
pC: ndarray[tuple[Any, ...], dtype[float64]]#
pE: ndarray[tuple[Any, ...], dtype[float64]]#
pF: ndarray[tuple[Any, ...], dtype[float64]]#
pY: ndarray[tuple[Any, ...], dtype[float64]]#
problem_type: str#
qp2lcp() → Tuple[LinearComplementarityProblem, Callable, Dict[str, Tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]]]][source]#

Transform LP/QP into LCP formulation, or transform pLP/pQP into pLCP formulation.

The original QP is in the following form:

(1)#\[\begin{split}\begin{aligned} \min\limits_x & \quad \frac{1}{2}x^\top Hx + (F \theta + f)^\top x \\ \text{s.t.} & \quad Ax \leq b + B\theta \\ & \quad A_e x = b_e + E\theta \end{aligned}\end{split}\]

There are 3 steps:

  1. Eliminate equality constraints in (1).

    To reduce the optimization problem to a simpler form, the system of linear equations \(A_e x = b_e + E\theta\) must be consistent, meaning no linearly dependent rows exist, and the number of equalities \(m_e\) is strictly less than the number of variables \(n\), i.e., \(m_e < n\).

    The principle is based on factorizing the equality constraints \(A_e x = b_e + E\theta\) into basic variables \(x_B\) and non-basic variables \(x_N\):

    \[A_{e,B} x_B + A_{e,N} x_N = b_e + E\theta,\]

    where the index sets \(B\) and \(N\) denote the columns for the factored system. The factored submatrix \(A_{e,B}\) must be invertible to express basic variables as a function of non-basic variables:

    \[x_{B} = - A_{e,B}^{-1} A_{e,N} x_N + A_{e,B}^{-1} b_e + A_{e,B}^{-1} E\theta.\]

    Using the substitutions:

    \[C = - A_{e,B}^{-1} A_{e,N}, D_1 = A_{e,B}^{-1} b_e, D_2 = A_{e,B}^{-1} E,\]

    the relation between basic and non-basic variables simplifies to:

    \[x_B = C x_N + D_1 + D_2 \theta.\]

    The original QP can be expressed with only non-basic variables \(x_N\) as follows:

    (2)#\[\begin{split}\begin{aligned} \min\limits_{x_N} &\quad \frac{1}{2}x_N \tilde{H} x_N + (\tilde{F} \theta + \tilde{f})^\top x_N \\ \text{s.t.} & \quad \tilde{A} x_N \leq \tilde{b} + \tilde{B}\theta \end{aligned}\end{split}\]

    with the following definitions (\(A_B\) and \(A_N\) are the columns of \(A\) corresponding to \(x_B\) and \(x_N\) respectively):

    \[\begin{split}\begin{aligned} & \tilde{H} = C^\top H_{BB} C + C^\top H_{BN} + H_{NB}C + H_{NN}, \\ & \tilde{f} = \frac{1}{2} C^\top \left(H_{BB}+H_{BB}^\top \right) D_1 + \frac{1}{2} \left(H_{BN}^\top + H_{NB}\right) D_1 + C^\top f_{B} + f_{N}, \\ & \tilde{F} = \frac{1}{2} C^\top \left(H_{BB}+H_{BB}^\top \right) D_2 + \frac{1}{2} \left(H_{BN}^\top + H_{NB}\right) D_2 + C^\top F_{B} + F_{N}, \\ & \tilde{A} = A_B C + A_N, \\ & \tilde{b} = b - A_B D_1, \\ & \tilde{B} = B - A_B D_2. \\ \end{aligned}\end{split}\]

  2. Transform the QP with only inequalities (2) into affine and non-negative form:

    (3)#\[\begin{split}\begin{aligned} \min\limits_y & \quad \frac{1}{2}y^\top \hat{H} y + (\hat{F} \theta + \hat{f})^\top y \\ \text{s.t.} & \quad \hat{A}y \geq \hat{b} + \hat{B}\theta \\ &\quad y \geq 0 \end{aligned}\end{split}\]

    Case 1: if \(\tilde{A}\) is not full row rank, factorize it row-wise into \(\tilde{A}_B\) and \(\tilde{A}_N\) where \(\tilde{A}_B\) is a square invertible matrix:

    \[\begin{split}\begin{aligned} -\tilde{A}_B x_N &= -\tilde{b}_B - \tilde{B}_B \theta + y \\ -\tilde{A}_N x_N &\geq -\tilde{b}_N - \tilde{B}_N \theta \\ y &\geq 0 \end{aligned}\end{split}\]

    Apply substitution \(x_N = \tilde{A}_B^{-1} \left( -y + \tilde{b}_B + \tilde{B}_B \theta \right)\) into (2), we can express the problem into the form of (3) with:

    \[\begin{split}\begin{aligned} & \hat{H} = \tilde{A}_B^{-\top} \tilde{H} \tilde{A}_B^{-1}, & & \hat{F} = - \hat{H} \tilde{B}_B - \tilde{A}_B^{-\top} \tilde{F}, & & \hat{f} = - \hat{H} \tilde{b}_B - \tilde{A}_B^{-\top} \tilde{f}, \\ & \hat{A} = \tilde{A}_N \tilde{A}_B^{-1}, & & \hat{b} = \tilde{A}_N \tilde{A}_B^{-1} \tilde{b}_B - \tilde{b}_N, & & \hat{B} = \tilde{A}_N \tilde{A}_B^{-1} \tilde{B}_B - \tilde{B}_N. \end{aligned}\end{split}\]

    Case 2: the \(x_N\) can be expressed as the difference of two positive vectors, i.e., \(x_N = x_N^+ - x_N^-\) with \(x_N^+\geq 0, ~x_N^-\geq 0\). By denoting \(y\) as their concatenation, we get the corresponding parts in (3):

    \[\begin{split}\begin{aligned} & \hat{H} = \begin{bmatrix} \tilde{H} & -\tilde{H} \\ \tilde{H} & \tilde{H} \end{bmatrix}, & & \hat{F} = \begin{bmatrix} \tilde{F} & -\tilde{F} \end{bmatrix}, & & \hat{f} = \begin{bmatrix} \tilde{f} & -\tilde{f} \end{bmatrix}, \\ & \hat{A} = \begin{bmatrix} -\tilde{A} & \tilde{A} \end{bmatrix}, & & \hat{b} = -\tilde{b}, & & \hat{B} = -\tilde{B}. \end{aligned}\end{split}\]

  3. Transform to LCP.

    By introducing slack variables \(s = \hat{A}y-\hat{b}-\hat{B}\theta\), the KKT condition of (3) is given by:

    \[\begin{split}\begin{aligned} & \hat{H}y + \hat{F}\theta + \hat{f} -\hat{A}^\top\lambda - \nu = 0 & \text{stationary point} \\ & \lambda^\top s = \nu^\top y = 0 & \text{complementarity} \\ & s, y\geq 0 & \text{primal feasibility} \\ & \lambda, \nu \geq 0 & \text{dual feasibility} \end{aligned}\end{split}\]

    The corresponding LCP is:

    \[\begin{split}\begin{bmatrix} v \\ s \end{bmatrix} - \begin{bmatrix} \hat{H} & -\hat{A}^\top \\ \hat{A} & 0 \end{bmatrix} \begin{bmatrix} y \\ \lambda \end{bmatrix} = \begin{bmatrix} \hat{F} \\ -\hat{B} \end{bmatrix} \theta + \begin{bmatrix} \hat{f} \\ -\hat{b} \end{bmatrix}\end{split}\]

The relationship between the solution of LCP and original QP is:

\[\begin{split}\begin{aligned} & x = P \begin{bmatrix} x_B \\ x_N \end{bmatrix} = P \left( \begin{bmatrix} C \\ I \end{bmatrix} x_N + \begin{bmatrix} D_1 \\ 0 \end{bmatrix} + \begin{bmatrix} D_2 \\ 0 \end{bmatrix}\theta \right), \\ & x_N = \begin{bmatrix} I & -I \end{bmatrix} y, \\ & y = \begin{bmatrix} I & 0 \end{bmatrix} z, \end{aligned}\end{split}\]

where \(P\) is a proper permutation matrix.

solve()[source]#

Solve the QP/pQP.

solve_for_param(theta: ndarray[tuple[Any, ...], dtype[float64]], **kwargs)[source]#

Solve the LCP / pLCP for a specific parameter value

ub: ndarray[tuple[Any, ...], dtype[float64]]#
mpt4py.solvers.qp.solve_pqp_with_ppopt(pqp: QuadraticProgram) → PqpSolution[source]#

Solve the parametric QP using PPOPT package.

In PPOPT, the standard form of quadratic multiparametric programming is defined as follows (from mpqp_program.py, https://ppopt.readthedocs.io/en/main/ppopt.html#module-ppopt.mpqp_program):

\[\begin{split}\begin{align} \min_x\quad \frac{1}{2}x^TQx& + \theta^TH^Tx + c^Tx\\ \text{s.t.}\quad Ax &\leq b + F\theta\\ A_{eq}x &= b_{eq}\\ A_\theta \theta &\leq b_\theta\\ x &\in R^n\\ \end{align}\end{split}\]

LCP and pLCP#

mpt4py.solvers.lcp.LCP#

alias of LinearComplementarityProblem

class mpt4py.solvers.lcp.LinearComplementarityProblem(M: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, q: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, Q: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, Ath: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>, bth: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.float64]] = <factory>)[source]#

Bases: object

Encapsulates data and solutions for LCP and pLCP.

\[w - Mz = q + Q\theta, w,z \geq 0, w^\top z = 0\]
Ath: ndarray[tuple[Any, ...], dtype[float64]]#
M: ndarray[tuple[Any, ...], dtype[float64]]#
Q: ndarray[tuple[Any, ...], dtype[float64]]#
bth: ndarray[tuple[Any, ...], dtype[float64]]#
d: int#
property is_parametric: bool#
n: int#
q: ndarray[tuple[Any, ...], dtype[float64]]#
solve(**kwargs) → LcpSolution | PlcpSolution[source]#

Solve the LCP / pLCP

solve_for_param(theta: ndarray[tuple[Any, ...], dtype[float64]], **kwargs) → LcpSolution[source]#

Solve the LCP / pLCP for a specific parameter value

LCP solver#

mpt4py.solvers.lcp_solve(M: ndarray[tuple[Any, ...], dtype[float64]], q: ndarray[tuple[Any, ...], dtype[float64]], solver: str = 'graves', max_pivots: int = 1000, max_time: float = 20, zero_tol: float = tolerances['zero'], verbose: bool = False) → LcpSolution[source]#

Solve a Linear Complementarity Problem (LCP)

pLCP solver#

mpt4py.solvers.plcp_solve(M: ndarray[tuple[Any, ...], dtype[float64]], Q: ndarray[tuple[Any, ...], dtype[float64]], q: ndarray[tuple[Any, ...], dtype[float64]], A_theta: ndarray[tuple[Any, ...], dtype[float64]] | None = None, b_theta: ndarray[tuple[Any, ...], dtype[float64]] | None = None, solver: str = 'lemke', zero_tol: float = tolerances['zero'], lex_tol: float = 1e-6, max_pivots: int = 1000, verbose: bool = False) → PlcpSolution[source]

Solve a Parametric Linear Complementarity Problem (pLCP).

Contents