Problem Formulation
PIQP expects QP of the form
\[\begin{aligned} \min_{x} \quad & \frac{1}{2} x^\top P x + c^\top x \\ \text {s.t.}\quad & Ax=b, \\ & Gx \leq h, \\ & x_{lb} \leq x \leq x_{ub} \end{aligned}\]with primal decision variables \(x \in \mathbb{R}^n\), matrices \(P\in \mathbb{S}_+^n\), \(A \in \mathbb{R}^{p \times n}\), \(G \in \mathbb{R}^{m \times n}\), and vectors \(c \in \mathbb{R}^n\), \(b \in \mathbb{R}^p\), \(h \in \mathbb{R}^m\), \(x_{lb} \in \mathbb{R}^n\), and \(x_{ub} \in \mathbb{R}^n\).
PIQP can handle infinite box constraints well, i.e. when elements of \(x_{lb}\) or \(x_{ub}\) are \(-\infty\) or \(\infty\), respectively. On the contrary, infinite values in the general inequalities \(Gx \leq h = \pm \infty\) can cause problems. PIQP internally disables them by setting the corresponding rows in \(G\) to zero (sparsity structure is preserved). For best performance, consider removing the corresponding constraints from the problem formulation directly.
Example QP
In the following the C++ interface of PIQP will be introduced using the following example QP problem:
\[\begin{aligned} \min_{x} \quad & \frac{1}{2} x^\top \begin{bmatrix} 6 & 0 \\ 0 & 4 \end{bmatrix} x + \begin{bmatrix} -1 \\ -4 \end{bmatrix}^\top x \\ \text {s.t.}\quad & \begin{bmatrix} 1 & -2 \end{bmatrix} x = 1, \\ & \begin{bmatrix} 1 & -1 \\ 2 & 0 \end{bmatrix} x \leq \begin{bmatrix} 0.2 \\ -1 \end{bmatrix}, \\ & -1 \leq x_1 \leq 1. \end{aligned}\]Problem Data
PIQP supports dense and sparse problem formulations. For small and dense problems the dense interface is preferred since vectorized instructions and cache locality can be exploited more efficiently, but for sparse problems the sparse interface and result in significant speedups.
The C++ interface uses Eigen and is automatically included in the PIQP header files
#include <limits>
#include "piqp/piqp.hpp"
<limits>
is included for getting infinity values.
We can then define the problem data as
int n = 2;
int p = 1;
int m = 2;
Eigen::MatrixXd P(n, n); P << 6, 0, 0, 4;
Eigen::VectorXd c(n); c << -1, -4;
Eigen::MatrixXd A(p, n); A << 1, -2;
Eigen::VectorXd b(p); b << 1;
Eigen::MatrixXd G(m, n); G << 1, -1, 2, 0;
Eigen::VectorXd h(m); h << 0.2, -1;
Eigen::VectorXd x_lb(n); x_lb << -1, -std::numeric_limits<double>::infinity();
Eigen::VectorXd x_ub(n); x_ub << 1, std::numeric_limits<double>::infinity();
For the sparse interface \(P\), \(A\), and \(G\) have to be in compressed sparse column (CSC) format.
Eigen::MatrixXd P_dense(n, n); P_dense << 6, 0, 0, 4;
Eigen::SparseMatrix<double> P = P_dense.sparseView();
Eigen::MatrixXd A_dense(p, n); A_dense << 1, -2;
Eigen::SparseMatrix<double> A = A_dense.sparseView();
Eigen::MatrixXd G_dense(m, n); G_dense << 1, -1, 2, 0;
Eigen::SparseMatrix<double> G = G_dense.sparseView();
There are better options to build sparse matrices in Eigen. For example, we can also build the sparse matrix \(P\) using the
insert
API:Eigen::SparseMatrix<double> P(n, n); P.insert(0, 0) = 6; P.insert(1, 1) = 4; P.makeCompressed();
For more details see the Eigen Docs.
Solver Instantiation
You can instantiate a solver object using
// for dense problems
piqp::DenseSolver<double> solver;
// or for sparse problems
piqp::SparseSolver<double> solver;
where the template argument defines the data type, i.e., double
or float
.
Settings
Settings can be directly set on the solver object:
solver.settings().verbose = true;
solver.settings().compute_timings = true;
In this example we enable the verbose output and internal timings. The full set of configuration options can be found here.
Solving the Problem
We can now set up the problem using
solver.setup(P, c, A, b, G, h, x_lb, x_ub);
Every variable except P
and c
are optional and laopt::nullopt
may be passed.
The data is internally copied, and the solver initializes all internal data structures.
Now, the problem can be solver using
piqp::Status status = solver.solve();
Status code
Status Code | Value | Description |
---|---|---|
PIQP_SOLVED | 1 | Solver solved problem up to given tolerance. |
PIQP_MAX_ITER_REACHED | -1 | Iteration limit was reached. |
PIQP_PRIMAL_INFEASIBLE | -2 | The problem is primal infeasible. |
PIQP_DUAL_INFEASIBLE | -3 | The problem is dual infeasible. |
PIQP_NUMERICS | -8 | Numerical error occurred during solving. |
PIQP_UNSOLVED | -9 | The problem is unsolved, i.e., solve was never called. |
PIQP_INVALID_SETTINGS | -10 | Invalid settings were provided to the solver. |
Extracting the Result
The result of the optimization can be obtained from the solver.result()
object. More specifically, the most important information includes
solver.result().x
: primal solutionsolver.result().y
: dual solution of equality constraintssolver.result().z
: dual solution of inequality constraintssolver.result().z_lb
: dual solution of lower bound box constraintssolver.result().z_ub
: dual solution of upper bound box constraintssolver.result().info.primal_obj
: primal objective valuesolver.result().info.run_time
: total runtime
Timing information like solver.result().info.run_time
is only measured if solver.settings().compute_timings
is set to true
.
Efficient Problem Updates
Instead of creating a new solver object everytime it’s possible to update the problem directly using
solver.update(P, c, A, b, G, h, x_lb, x_ub);
with a subsequent call to
piqp::Status status = solver.solve();
This allows the solver to internally reuse memory and factorizations speeding up subsequent solves. Similar to the setup
function, all parameters are optional and laopt::nullopt
may be passed instead.
Note the dimension and sparsity pattern of the problem are not allowed to change when calling the update
function.