Antoine Bambade
email: author.bambade.antoine@gmail.com

Fabian Schramm

Sarah El Kazdadi

Stéphane Caron

Adrien Taylor

Justin Carpentier

Sarah El-Kazdadi

PROXQP: an Efficient and Versatile Quadratic Programming Solver for Real-Time Robotics Applications and Beyond

Keywords: Optimization, Quadratic Programming, Embedded optimization

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Optimization has become a key enabler to simplify and systematize the programming of complex robot motions. Nowadays, many robotic problems, ranging from simulation and control to planning and estimation, are framed as optimization problems. An important class of such optimization problems is that of convex quadratic programs (QP), which allows dealing with friction-less unilateral contact modelling [START_REF] Redon | Gauss least constraints principle and rigid body simulations[END_REF], constrained forward dynamics [START_REF] Carpentier | Proximal and sparse resolution of constrained dynamic equations[END_REF], inverse kinematics and dynamics for task control [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF]- [START_REF] Kuindersma | Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot[END_REF], and legged locomotion [START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF]- [START_REF] Wensing | Optimization-based control for dynamic legged robots[END_REF], among others. QPs are also commonly used as subroutines for solving more complex problems, for instance, in the context of constrained optimal control problems [START_REF] Leineweber | An efficient multiple shooting based reduced SQP strategy for large-scale dynamic process optimization. part 1: theoretical aspects[END_REF]- [START_REF] Jallet | Constrained differential dynamic programming: A primal-dual augmented lagrangian approach[END_REF].

Formally, a QP corresponds to the minimization of a convex quadratic cost under linear inequality constraints. It is mathematically described as follows:

min x∈R n f (x) def = 1 2 x ⊤ Hx + g ⊤ x s.t. Cx ⩽ u, (QP)
where H ∈ R n×n is a real symmetric positive semi-definite matrix (notation S n + ), g ∈ R n , C ∈ R m×n , and u ∈ R m . The problem dimension is n, while m corresponds to the numbers inequality constraints.

In many scenarios, QP instances have to be solved at high frequency (e.g., 1 kHz is typical for inverse kinematics or dynamics), under various levels of accuracy depending on the application, and potentially in relatively large dimensions (e.g., humanoid robots) particularly when it comes to model predictive control (MPC). In such contexts, it is valuable to exploit as much as possible the structure offered by the task (e.g., sparsity pattern of the underlying problem) and to make use of an optimization method that benefits from various warmstarting features. Yet, in a robotic context, warm-starting a QP with an initial guess from a previous instance requires it to lie in the feasible set of the new instance. This property, also called recursive feasibility [START_REF] Ciocca | Strong recursive feasibility in model predictive control of biped walking[END_REF] in the MPC literature, is challenging to enforce for closed-loop systems as the latest state estimated from sensory measurements may not be feasible, prompting practitioners to relax the initial state constraint [START_REF] Mayne | Robust model predictive control of constrained linear systems with bounded disturbances[END_REF], [START_REF] Villa | Model predictive control of biped walking with bounded uncertainties[END_REF]. These practical reasons may limit the use of warm-starting and hotstarting to sequential problems for which recursive feasibility can be guaranteed, such as open-loop [START_REF] Ciocca | Strong recursive feasibility in model predictive control of biped walking[END_REF] rather than closedloop [START_REF] Villa | Model predictive control of biped walking with bounded uncertainties[END_REF] MPC.

This work develops a new augmented Lagrangian-based QP method and solver with the following contributions: → The PROXQP algorithm itself, a generic method for solving convex QPs, based on a primal-dual proximal augmented Lagrangian method. In particular, we provide a convergence proof of the overall algorithm and highlight that the PROXQP algorithm, and more generally primaldual Proximal augmented Lagrangian methods, actually solves the closest primal-feasible QP-in a classical ℓ 2 sense that we detail in the sequel-if the original QP appears to be primally infeasible. → On the software side, we distribute an open-source and flexible implementation of PROXQP within the PROXSUITE library https://github.com/Simple-Robotics/proxsuite. The solver is written in C++ for efficiency and comes with Julia and Python bindings. → On the numerical side, we illustrate the practical performance of PROXQP on different standard sparse and dense robotics and control experiments, including realworld closed-loop problems. We notably highlight how solving the closest primal feasible QPs can be leveraged in closed-loop model predictive control. We also benchmark the performance of PROXQP against state-of-the-art QP solvers on various problem sets of generic QPs from the literature.

This article is an extension of the previously published conference paper [START_REF] Bambade | PROX-QP: Yet another Quadratic Programming Solver for Robotics and beyond[END_REF]. The revised version of PROXQP we derive here is more efficient and handles sparse problems. Additionally, we now prove that the PROXQP algorithm features a global convergence guarantee, as well as other advantageous numerical properties (see Section III-B). Moreover, we highlight a central feature of PROXQP for adequately handling primally infeasible QPs, which can be exploited in a closed-loop model predictive control scheme, as illustrated in the experimental section (see Sec. VI) or in machine learning context as leveraged in [START_REF] Bambade | QPLayer: efficient differentiation of convex quadratic optimization[END_REF]. Finally, this article extends the validation benchmark to include standard control problems and a variety of robotics tasks.

The remainder of this paper is organized as follows. Section II reviews augmented Lagrangian methods for solving quadratic programs, which are at the core of our approach. Section III presents the PROXQP algorithm and its theoretical features. Section IV provides an overview of the diverse optimization methods commonly used for solving QPs, and highlights the benefits of the primal-dual augmented Lagrangian approach leveraged in PROXQP. Section V details the various features of the PROXQP implementation within the ProxSuite library. Finally, Section VI provides an extensive validation of PROXQP on different robotic and control tasks. We also benchmark the numerical robustness, accuracy, and speed of our approach against state-of-the-art solvers on more generic QPs. Due to space limitations and to facilitate the overall reading of the paper, the annexes are collected in the companion report [START_REF] Bambade | Companion Report of PROXQP: an Efficient and Versatile Quadratic Programming Solver for Real-Time Robotics Applications and Beyond[END_REF].

II. BACKGROUND

In this section, we recall some of the well-established components underlying PROXQP, namely, (i) optimality conditions related to (QP), and (ii) existing augmented Lagrangian-based strategies. The reader already familiar with this topic can skip this section and directly dive into Section III, which details the PROXQP algorithm.

A. Notation

The Lagrangian L associated to (QP) is defined by: L(x, z)

def =f (x) + z ⊤ (Cx -u),
where x ∈ R n are the primal variables and z ∈ R m + corresponds to the dual variables. The dual function is then defined for z ⩾ 0 as δ(z) def = inf x∈R n L(x, z). The dual problem to (QP) corresponds to: max z⩾0 δ(z).

(Dual-QP)

We say that (Dual-QP) is feasible when δ < +∞.

B. Optimality conditions

For linearly constrained convex optimization problems such as (QP), strong duality holds (when either the primal or the dual problem is feasible), and the associated Karush-Kuhn-Tucker (KKT) conditions are necessary and sufficient for optimality of a primal-dual point (x, z) (see e.g. [START_REF] O'donoghue | Conic optimization via operator splitting and homogeneous self-dual embedding[END_REF]Section 2.1]). For (QP), the KKT system is:

       Hx + g + C ⊤ z = 0, z ⩾ 0, Cx ⩽ u, x ⊤ Hx + g ⊤ x + u ⊤ z = 0. (KKT)
In practice, we look for a tuple (x, z) satisfying these optimality conditions (KKT) up to a certain level of predefined accuracy ϵ abs > 0 (whose value typically depends on the application at hand), in the following sense:

   ∥Hx + g + C ⊤ z∥ ∞ ⩽ ϵ abs , ∥[Cx -u] + ∥ ∞ ⩽ ϵ abs , |x ⊤ Hx + g ⊤ x + u ⊤ z| ⩽ ϵ abs , (1)
where [.] + stands for the component-wise projection onto the nonnegative orthant, i.e. [v] +,i = max(0, v i ). In this work, we prefer the ℓ ∞ norm to the ℓ 2 norm as it does not scale with the dimension d. It is also common practice to rely on relative error criteria for early-stopping, as absolute targets might be out of reach due to numerical issues [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]- [START_REF] Optimization | GUROBI optimizer reference manual[END_REF]. Although this option is included in the software, in what follows, we will not add relative terms after absolute ones for the sake of simplicity.

C. Method of Multipliers

The method of multipliers (MM) is based on the minimization of the augmented Lagrangian (AL) penalty function, and is a precursor to the algorithm developed in this work. ALbased approaches typically feature advantageous properties such as convergence toward a solution to the closest feasible QP (see Section VI for examples in robotics). Augmented Lagrangian penalty function. The so-called AL function [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Section 4] relies on augmenting the cost function f with a shifted ℓ 2 penalization of the constraint (i.e., Cx -u) with respect to a weighted dual multipliers µz:

LA(x, z; µ) def = f (x) + 1 2µ ∥[Cx -u + µz]+∥ 2 2 -∥µz∥ 2 2 ,
where µ > 0 is a positive penalty parameter. In the case where we replace the inequality constraints Cx ⩽ u of (QP) by equalities Cx = u, this would be equivalent to augmenting the standard Lagrangian L with a weighted penalization:

L A (x, z; µ) = L(x, z) + 1 2µ ∥Cx -u∥ 2 2 .
Method of multipliers. Many optimization strategies rely on alternating (i) minimizing L A with respect to the primal variables, and (ii) maximizing L A with respect to the dual variables. MM [START_REF] Hestenes | Multiplier and gradient methods[END_REF], [START_REF] Powell | A method for nonlinear constraints in minimization problems[END_REF] is one such strategy, where the maximization step (ii) is computed in closed-form. More precisely, considering some sequence {µ k } k bounded below by some µ ∞ > 0:

x k+1 = arg min x∈R n L A (x, z k ; µ k ), z k+1 = arg max z∈R m L A (x k+1 , z; µ k ) = z k + 1 µ k (Cx k+1 -u) + .
Yet, in the presence of inequality constraints, the exact minimization of L A (x, z k ; µ k ) with respect to x is generally hard to compute in closed-form since x → L A (x, z k ; µ k ) is piece-wise quadratic [START_REF] Sun | On piecewise quadratic Newton and trust region problems[END_REF]. For this reason, practical MM solvers typically rely on ϵ k -approximate solutions for computing x k+1 :

x k+1 ≈ ϵ k arg min x∈R n L A (x, z k ; µ k ), z k+1 = z k + 1 µ k (Cx k+1 -u) + , (2)
for some notion of "≈ ϵ k ". For instance, we can use:

∥∇ x L A (x k+1 , z k ; µ k )∥ ∞ ⩽ ϵ k .
For ensuring convergence of the numerical scheme, one typically needs to enforce certain properties on ϵ k and {µ k } k .

For instance, it is clear that one must have ϵ k → 0; other conditions include summability of the sequence {ϵ k } k , see, e.g., [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF], [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF], [START_REF] Nocedal | Numerical optimization[END_REF]. The boundedness condition on {µ k } k turns out also to be necessary for ensuring convergence [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF].

Convergence properties. It turns out that MM can be seen as a proximal point method applied to (Dual-QP), see [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Section 4]. Therefore, it is guaranteed that {z k } k converges to an optimal solution to (QP) as soon as there exists one [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Theorem 4]. Furthermore, {x k } k is a minimizing sequence, i.e., f (x k ) -→ min x∈R n f (x). Note that {x k } k does not necessarily converge towards a solution to (QP) [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]. This equivalence between MM and the dual proximal point algorithm also allows MM to inherit a key property for handling primal infeasible problems. Indeed, assuming (Dual-QP) to be feasible, the dual iterates of MM actually converge linearly to a solution of the dual to the following hierarchical problem [START_REF] Chiche | How the augmented Lagrangian algorithm can deal with an infeasible convex quadratic optimization problem[END_REF]:

s ⋆ = arg min s∈R m 1 2 ∥s∥ 2 2 s.t. x ⋆ , z ⋆ ∈ arg min x∈R n max z∈R m + L(x, z, s), (QP-H) with L(x, z, s) def = f (x) + z ⊤ (Cx -u -s)
(namely the Lagrangian of (QP) augmented with a shift variable s). Under strict convexity of f , it results that as soon as (QP) is primal infeasible (i.e., there is no x s.t. Cx ⩽ u), then MM will converge towards a solution (x ⋆ , z ⋆ ) of (QP-H). Such a solution satisfies the following KKT optimality conditions:

Hx ⋆ + g + C ⊤ z ⋆ = 0, Cx ⋆ ⩽ u + s ⋆ , (x ⋆ ) ⊤ Hx ⋆ + g ⊤ x ⋆ + (u + s ⋆ ) ⊤ z ⋆ = 0.
In other words, s is the smallest possible (measured in ℓ 2 norm) shift on the constraints that makes (QP) feasible.

D. Proximal Method of Multipliers

The proximal method of multipliers (PMM) is an extension of MM with stronger convergence guarantees. Proximal augmented Lagrangian function. For a certain ρ > 0 ( 1 /ρ corresponds to a step-size), PMM relies on an additional proximal term on the primal variables with [23, Equation 1.9], considering an additional sequence {µ k } k bounded below by some µ ∞ > 0:

(x k+1 , z k+1 ) = arg min x∈R n max z∈R m + L(x, z) + ρ 2 ∥x -x k ∥ 2 2 -µ k 2 ∥z -z k ∥ 2 2 .
(

) 3
Assuming that one can solve (3), one can generate a sequence {(x k+1 , z k+1 )} k that converges to a primal-dual optimal pair for (QP). However, (3) involves a convex-concave saddle-point problem, for which numerous numerical methods exist yet closed-form solutions can seldom be found. For solving this problem, a natural strategy consists in splitting the optimization procedure in two parts:

x k+1 = arg min x∈R n Φ k ρ,µ k (x), z k+1 = arg max z∈R m + L(x, z) -µ k 2 ∥z -z k ∥ 2 2 ,
where

Φ k ρ,µ k (x) def = L A (x, z k ; µ k )+ ρ 2 ∥x-x k ∥ 2
2 is often referred to as the Proximal augmented Lagrangian [START_REF] Hermans | Qpalm: A newton-type proximal augmented lagrangian method for quadratic programs[END_REF], and where the second line can be computed in closed form as

z k+1 = z k + 1 µ k (Cx k+1 -u)] + .
As with MM, one generally needs to resort on numerical methods for the first minimization stage, which we can only solve approximately:

x k+1 ≈ ϵ k arg min x∈R n Φ k ρ,µ k (x), z k+1 = z k + 1 µ k (Cx k+1 -u) + ,
for some notion of "≈ ϵ k ". For instance, we can use:

∇ x Φ k ρ,µ k (x k+1 ) ∞ ⩽ ϵ k .
Again, it is clear that convergence of {(x k , z k )} towards a primal-dual solution requires certain properties of the approximation error, such as ϵ k → 0. Other conditions include summability of {ϵ k } k , and {µ k } being bounded from below by some µ ∞ > 0, see, e.g., [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Theorem 5].

Convergence properties. PMM enjoys the key advantage of guaranteed primal-dual convergence of its iterates under relatively weak assumptions, namely the existence of an optimal primal-dual pair with zero duality gap [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Theorem 7]. This guarantee is stronger than that of MM, for which only the dual variables are guaranteed to converge [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Theorem 4]. Furthermore, the proximal augmented Lagrangian function is a piece-wise quadratic strongly convex function, unlike the augmented Lagrangian x → L A (x, z k ; µ k ) which is only guaranteed to be convex. This allows accelerating the computation of the iterates, e.g., using quadratic Newton-type methods [START_REF]QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs[END_REF]. We leverage this key feature in PROXQP, as further detailed in Section III. Finally, PMM inherits from MM convergence towards the closest feasible QP solution when the problem is primal infeasible (QP-H), see Remark 1.

E. A primal-dual Proximal Method of Multipliers

Primal-dual proximal methods of multipliers (PDPMM) is an alternative formulation of PMM. Its name stems from its choice of a primal-dual penalty function for evaluating the quality of the iterates of the sequence of proximal sub-problems. When intermediary optimization subproblems are solved exactly, PMM and PDPMM produce the same iterates. Thereby, PDPMM benefits from the strong convergence properties of PMM. However, the alternate formulation of PMM is convenient and can be leveraged at computation time, and we leverage it in our algorithms. Primal-dual proximal augmented Lagrangian function. The author of [START_REF] Marchi | On a primal-dual newton proximal method for convex quadratic programs[END_REF] shows that (3) can be equivalently reformulated as a minimization of the following primal-dual merit function:

M k ρ,µ (x, z) def =f (x) + 1 µ [Cx -u + µ(z k -z 2 )] + 2 2 + ρ 2 ∥x -x k ∥ 2 2 + µ 4 ∥z∥ 2 2 . (4)
Hence ( 3) is equivalent to:

(x k+1 , z k+1 ) = arg min x∈R n ,z∈R m M k ρ,µ k (x, z). (5)
Of course, solving this optimization step still requires using numerical schemes, and can therefore only be approximated:

(x k+1 , z k+1 ) ≈ ϵ k arg min x∈R n ,z∈R m M k ρ,µ k (x, z), (6)
for some notion of "≈ ϵ k ". For instance, [START_REF] Marchi | On a primal-dual newton proximal method for convex quadratic programs[END_REF] uses:

∇ x f (x k+1 ) + C ⊤ z k+1 + ρ(x k+1 -x k ) [Cx k+1 -u + µ k (z k -z k+1 /2)] + -µ k z k+1 /2 ∞ ⩽ ϵ k .
The PROXQP algorithm relies on a similar error criterion as detailed in the next sections. Convergence properties. As for PMM, the recurrence (6) guarantees that {(x k , z k )} k converges to an optimal primaldual pair under the relatively weak assumptions (the same as in (II-D0a): as soon as there exists a solution to (QP), and {µ k } k is bounded below by some µ ∞ > 0). Furthermore, M k ρ,µ k is both strongly convex and piece-wise quadratic with respect to both the primal and the dual variables. Hence, as for Φ k ρ,µ , its minimization can be performed efficiently using semi-smooth Newton-type methods. As we show in Section III, the linear systems induced by the minimization of M k ρ,µ k are naturally better conditioned than those arising from Φ k ρ,µ k , which is due to the fact that this formulation naturally avoids quadratic matrix multiplications (such as C ⊤ C). Finally, Remark 1 highlights that more generally to MM methods, PMM and PDPMM both converge towards the primal-dual pair of the closest feasible QP solution, which is useful when (QP) is infeasible.

Remark 1 (Convergence towards solutions of (QP-H)). It can be proven following [START_REF] Chiche | How the augmented Lagrangian algorithm can deal with an infeasible convex quadratic optimization problem[END_REF] that if (QP) is primal infeasible, its dual (Dual-QP) is feasible, and {µ k } k is a sequence satisfying ∀k ∈ N, µ k ⩾ µ ∞ for some µ ∞ > 0, then, the sequence {(x k , z k )} generated by (3) converge towards a solution of (QP-H).

F. Discussions

There are three practical ingredients that are key when implementing AL-based approaches. 1 -Setting ϵ k . Appropriate choices for ϵ k are fundamental to any AL method.

• On the one hand, picking too small values for ϵ k leads to unnecessary good approximations of the intermediate optimization problems, even when the iterates are still far away from solutions to (QP). Practically, this scenario corresponds to requiring the internal solver to unnecessarily run for too long when solving the intermediary problems. • On the contrary, picking ϵ k too large results in faster computations of the intermediary subproblems at the cost of requiring much more iterations of the AL method, which might even possibly not converge.

There are many different techniques for fixing ϵ k in practice [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF], [START_REF]QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs[END_REF], [START_REF] Birgin | Practical Augmented Lagrangian Methods for Constrained Optimization[END_REF]. For PROXQP, we chose to rely on an adaptive strategy similar in spirit to [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF].

2 -Setting µ k . Using different values for µ k has a direct impact on the practical convergence speed (along with the value ρ for PMM and PDPMM) (see, e.g., [START_REF] Aspremont | Acceleration methods[END_REF]Chapter 5])

• On the one hand, lower values of µ k lead to faster convergence (i.e., fewer iterations). • On the other hand, larger values of µ k render the intermediary subproblems structurally simpler to solve (by increased strong convexity).

Practically speaking, let us mention that in many solvers (including PROXQP), updating µ k corresponds to re-factorizing some internal matrices. Hence, we must generally limit the number of updates to µ k for best practical performances.

3 -Efficient method for computing (x k+1 , z k+1 ). There are several different ways to come up with ϵ k -approximations for (x k+1 , z k+1 ). The choice of the internal routine for this approximation is thereby crucial both for timing and accuracy requirements.

The next section presents our AL-based approach for solving (QP), which we refer to as PROXQP. In particular, we discuss how to handle these three key aspects for computational efficiency.

III. PROXQP

In this section, we introduce the PROXQP algorithm for solving generic QP problems. For clarity issues, we provide the proofs of convergence of PROXQP in the companion report [START_REF] Bambade | Companion Report of PROXQP: an Efficient and Versatile Quadratic Programming Solver for Real-Time Robotics Applications and Beyond[END_REF].

A. The PROXQP algorithm PROXQP is detailed in Algorithm 1. It combines a proximal augmented Lagrangian technique (approximated using an internal primal-dual method) with a globalization strategy for scheduling the internal step-size and accuracy parameters. Inputs:

• initial states:

x 0 , z 0 , • initial parameters: ϵ 0 bcl , ϵ 0 , ϵ abs , ρ, µ 0 > 0 • hyper-parameters: µ f , β bcl ∈ (0, 1), α bcl ∈ (0, 1/2) with α bcl + β bcl < 1, µ min > 0, k max ∈ N ∪ {+∞}. Initialization: • preconditioning (see paragraph V-B0d) • optional initialization (see paragraph V-B0b) of x 0 ,z 0 . while Stopping criterion (1) not satisfied do
Compute (x k+1 , ẑk+1 ) satisfying ( 11) at accuracy ϵ k (using a semi-smooth Newton method, see Section III-D);

if p k ⩽ ϵ k bcl using (8) OR k ⩾ k max then µ k+1 = µ k , ϵ k+1 = ϵ k µ k+1 , ϵ k+1 bcl = ϵ k bcl (µ k+1 ) βbcl z k+1 = ẑk+1 else µ k+1 = max(µ min , µ f µ k ) ϵ k+1 = ϵ 0 µ k+1 , ϵ k+1 bcl = ϵ 0 bcl (µ k+1 ) α bcl z k+1 = z k end k ← k + 1 end Output: A (x k , z k ) satisfying the ϵ abs -approximation criterion (1) for problem (QP).
More precisely, at iteration k, PROXQP computes candidates x k+1 , ẑk+1 satisfying

(x k+1 , ẑk+1 ) ≈ ϵ k arg min x∈R n ,z∈R m M k ρ,µ k ,α (x, z), (7)
for a merit function M k ρ,µ k ,α that we detail in Section III-D. The following steps build on the bound-constrained Lagrangian (BCL) strategy (see, e.g., [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF] and [START_REF] Nocedal | Numerical optimization[END_REF]Algorithm 17.4]) for scheduling both ϵ k and µ k .

The main idea consists in a fast decrease of ϵ k when the primal residual ∥[C xk+1 -u] + ∥ ∞ is small enough. In this case, we accept the approximation with z k+1 = ẑk+1 and leave µ k unchanged with µ k+1 = µ k . Otherwise, we slightly decrease ϵ k , reject the approximation with z k+1 = z k and decrease µ k . Decreasing µ k corresponds to a heavier penalization of the feasibility constraint within the augmented Lagrangian-based function at the next step. Such a strategy has been proved to perform well in advanced optimization packages such as LANCELOT [START_REF] Conn | Lancelot a Fortran Package for Large-Scale Nonlinear Optimization[END_REF] as well as in robotics for solving constrained optimal control problems [START_REF] Jallet | Constrained differential dynamic programming: A primal-dual augmented lagrangian approach[END_REF], [START_REF] Plancher | Constrained unscented dynamic programming[END_REF]- [START_REF] Jallet | Implicit differential dynamic programming[END_REF]. We detail this globalization technique alongside its resulting global convergence properties for PROXQP in Section III-B. The method used to minimize M k ρ,µ k ,α is explained in Section III-D.

B. Globalization strategy for scheduling ϵ k and µ k

The globalization strategy in PROXQP consists in updating the dual variables z k -obtained from minimizing a primaldual proximal augmented Lagrangian merit function (described Algorithm 2: PROXQP (analyzed version)

Inputs:

• initial states:

x 0 , z 0 , • initial parameters: ϵ 0 bcl , ϵ 0 > 0, ρ 0 = µ 0 ∈ (0, 1) • hyper-parameters: µ f < 1, α bcl ∈ (0, 1/2),β bcl ∈ (0, 1) with α bcl + β bcl < 1, k max ∈ N ∪ {+∞}. for k = 0, 1, . . . do Compute (x k+1 , ẑk+1 ) satisfying (11) at accuracy ϵ k with ρ k = µ k ; if p k ⩽ ϵ k bcl OR k ⩾ k max then µ k+1 = µ k , ϵ k+1 = ϵ k µ k+1 , ϵ k+1 bcl = ϵ k bcl (µ k+1 ) βbcl , z k+1 = ẑk+1 else µ k+1 = µ f µ k , ϵ k+1 = ϵ 0 µ k+1 , ϵ k+1 bcl = ϵ 0 bcl (µ k+1 ) αbcl , z k+1 = z k . end end
in Section III-D)-only when the corresponding primal inequality constraint violation (denoted by p k hereafter) is small enough. More precisely, we use a second internal sequence of tolerances denoted by ϵ k bcl (which is also tuned within this BCL strategy) and update the dual multipliers only when p k ⩽ ϵ k bcl , where p k denotes the ℓ ∞ primal feasibility violation:

p k def = [Cx k+1 -u] + ∞ . (8)
This BCL strategy allows searching for appropriate values for the hyper-parameters ϵ k , ϵ k bcl , µ k . For the constraint penalization parameter µ k , it proceeds as follows:

• If p k ⩽ ϵ k
bcl : the primal inequality constraint violation is satisfactory, keep µ k as is.

• Otherwise: the primal inequality constraint violation is too large, decrease µ k for subsequent proximal subproblems, which has the effect of leading to heavier penalization of infeasibility, enforcing their satisfaction.

About the accuracy parameters ϵ k and ϵ k bcl , the update rules are more technical. The motivation underlying those choices is to ensure global convergence: a geometric-decay type update when the primal inequality constraint violation is good enough, and see [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF]Lemma 4.1] for when the infeasibility is too large. The detailed strategy is summarized in Algorithm 1. Convergence properties. We analyze a simpler algorithm (see Algorithm 2), whose small variations are more effective in practice. More precisely, it assumes that at each iteration k of the algorithm ρ k = µ k compared to Algorithm 1, where ρ is fixed. Updating ρ is not efficient in practice since it would require a numerical factorization at each iteration. In this setup and under the assumption that there exists a primaldual solution to (QP) (see Theorem 1 and Remark 2), the iterates of PROXQP are guaranteed to converge to an optimal primal-dual pair.

Assumption 1. The problem (QP) is feasible. Assumption 2. The iterates {(x k , z k )} k generated by PROXQP are bounded.
Theorem 1 (Convergence of PROXQP). Under Assumption 1 (existence of a solution to QP), Assumption 2 (bounded iterates), and with all parameters set as required in the inputs of Algorithm 2, the iterates {(x k , z k )} of PROXQP converge to a solution (x ⋆ , z ⋆ ) of (QP).

Proof. See (Appendix-E) Lemma 1. Under Assumption 1 (existence of solution to (QP)) and Assumption 2 (bounded iterates), and with all parameters set as required in the inputs of Algorithm 2, ∃K ∈ N, and

µ ∈ R such that ∀k ⩾ K µ k = µ. Hence, ∀k ⩾ 0, p K+k ⩽ ϵ K bcl (µ K ) β bcl k .
Remark 2 (On the necessity of Assumption 2). The BCL strategy [27, Algorithm 1 and Algorithm 2] was originally developed for solving more general optimization problems. BCL relies on an augmented Lagrangian strategy. In other words, we build on BCL for constructing an algorithm that exploits the structural properties of simpler optimization problems (namely QPs). Furthermore, rather than being based on a purely AL strategy, we rely on a proximal AL strategy. We managed to get rid of most of the strong assumptions required for ensuring global convergence of BCL [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF]. Yet, Assumption (2) still appeared necessary in our analysis. We introduced the safeguard parameter k max ∈ N to always enforce the algorithm to ultimately always accept the candidates (7) if Assumption 2 does not hold. This simple trick allows inheriting some nice properties from PMM that include convergence under mild assumptions [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Theorem 5].

Remark 3 (Advantageous numerical properties).

As a byproduct of the proof of Theorem 1, the method is guaranteed to use constant penalization parameters after a finite number of iterations. This is advantageous insofar as (1) we inherit favorable properties of proximal point methods on saddle point, and (2) it limits the number of numerical matrix factorizations.

C. Infeasible QPs

Finally, Corollary 1 guarantees that, similarly to MM, PMM and PDPMM methods, PROXQP converges towards the primaldual pair to the closest feasible QP solution, which is useful when (QP) is infeasible.

Corollary 1. Let (QP) be primal infeasible, let its dual (Dual-QP) be feasible, let µ > 0 and k max ∈ N. The sequence {(x k , z k )} generated by PROXQP with exact subproblem minimization converges towards a solution of (QP-H).

Proof. This is a direct consequence of the fact that, after a finite number of iterations, the iterations of PROXQP correspond to those of PMM.

D. Solving proximal sub-problems

We now explicit the primal-dual proximal augmented Lagrangian merit function used for generating candi-dates (x k+1 , ẑk+1 ) in [START_REF] Carpentier | Recent progress in legged robots locomotion control[END_REF]. Then, we explain how to minimize it in quadratic time. The PROXQP primal-dual proximal augmented Lagrangian merit function. We use the same procedure as in [START_REF] Marchi | On a primal-dual newton proximal method for convex quadratic programs[END_REF] for building our primal-dual proximal augmented Lagrangian merit function M k ρ,µ k ,α . By construction, it guarantees that minimizing M k ρ,µ k ,α is equivalent to solving (3). Yet, contrary to M ρ,µ k , we make a different choice for the parameter α ∈ (0, 1) 1 , tuned for better practical performance,

M k ρ,µ k ,α (x, z) def =f (x) + 1 2αµ k ∥[Cx -u + µ k (z k + (α -1)z)]+∥ 2 2 + ρ 2 ∥x -x k ∥ 2 2 + (1-α)µ k 2 ∥z∥ 2 2 . ( 9
)
This leads to the following inexact proximal scheme

(x k+1 , z k+1 ) ≈ ϵ k arg min x∈R n ,z∈R m M k ρ,µ k ,α (x, z), (10)
where

≈ ϵ k corresponds to ∇xf (x k+1 ) + C ⊤ z k+1 + ρ(x k+1 -x k ) [Cx k+1 -u + µ k (z k + (α -1)z k+1 )]+ -αµ k z k+1 ∞ ⩽ ϵ k . ( 11
)
Since M k ρ,µ k ,α is strongly convex and piece-wise quadratic, it can be minimized exactly in finite time using a semi-smooth Newton method with exact linesearch [26, Theorem 3]. Semi-smooth Newton method. For convenience, we omit the index k for the iterates of the Newton method arising at iteration k of PROXQP. A semi-smooth Newton method applied to M k ρ,µ k ,α , and initialized at (x (0) , ẑ(0) ) = (x k , z k ) involves finding at the l th inner iteration dw def = (dx, dz) satisfying:

∇ 2 M k ρ,µ k ,α (x (l) , ẑ(l) )dw + ∇M k ρ,µ k ,α (x (l) , ẑ(l) ) = 0, (12)
where ∇M k ρ,µ k ,α (x (l) , ẑ(l) ) and ∇ 2 M k ρ,µ k ,α (x (l) , ẑ(l) ) stand respectively for one element of the generalized Jacobian and the generalized Hessian of M ρ,µ k ,α [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF]Section 7.1].

Once a semi-smooth Newton step (dx, dz) has been obtained, the exact line-search procedure consists in finding the unique t * such that:

t * = arg min t⩾0 M k ρ,µ k ,α (x (l) + tdx, ẑ(l) + tdz).
Similarly to [START_REF] Bambade | PROX-QP: Yet another Quadratic Programming Solver for Robotics and beyond[END_REF], [START_REF]QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs[END_REF], [START_REF] Marchi | On a primal-dual newton proximal method for convex quadratic programs[END_REF] the function t → M k ρ,µ k ,α (x (l) + tdx, ẑ(l) + tdz) is continuous and piece-wise quadratic with a finite number of breaking points 2 . Hence, one can exactly compute t * , an optimal solution to this one-dimensional problem. Finally, the semi-smooth Newton method initialized at (x (0) , ẑ(0) ) = (x k , z k ) generates a sequence (x (l) , ẑ(l) ) via the update rule:

x(l+1) = x(l) + t * dx, ẑ(l+1) = ẑ(l) + t * dz. 1 Notice that with α = 1 2 we recover M ρ,µ k as in (4). 2 This feature is very specific to linear and quadratic types of problems and does help reducing the overall computational burden.

Solving linear systems. In this section, we explicit the linear systems induced by the minimization of M k ρ,µ k ,α . We show that they are naturally better conditioned than those arising from Φ k ρ,µ k since they naturally avoid quadratic matrix multiplications (such as C ⊤ C).

Equation [START_REF] Jallet | Constrained differential dynamic programming: A primal-dual augmented lagrangian approach[END_REF] reads:

H + ρI + 1 αµ k C ⊤ (I -P (k,l) )C α-1 α C ⊤ (I -P (k,l) ) α-1 α (I -P (k,l) )C µ k (α -1)[-1 α I -α-1 α P (k,l) ] dx dz = - ∇xf (x (l) ) + ρ(x (l) -x k ) + 1 αµ k C ⊤ [w (l) k -u] + α-1 α [w (l) k -u] + + µ k (1 -α)ẑ (l) , (13)
where P (k,l) stands for one element of the generalized Jacobian of [.] + at w l) . P (k,l) consists of a diagonal matrix with entries

(l) k def = Cx (l) + µ k z k + µ k (α -1)z (
P (k,l) jj def = 1 if C ⊤ j x(l) + µ k z k j + µ k (α -1)ẑ (l) j ⩽ u j , 0 otherwise. When P (k,l) jj = 1, it implies (see [32, Section 3.2] for more details) dz j = -(ẑ (l) ) j .
The linear system (13) can hence be equivalently formulated as:

               H + ρI C ⊤ J k (x (l) ,ẑ (l) ) C J k (x (l) ,ẑ (l) ) -µ k I J k (x (l) ,ẑ (l) ) dx dz J k (x (l) ,ẑ (l) ) = - H x(l) + g + ρ(x (l) -x k ) + C ⊤ J k (x (l) ,ẑ (l) ) ẑ(l) J k (x (l) ,ẑ (l) ) [Cx (l) + µ k z k -µ k ẑ(l) -u] J k (x (l) ,ẑ (l) ) , dz J k (x (l) ,ẑ (l) ) c = -ẑ (l) J k (x (l) ,ẑ (l) ) c , (14)
where:

J k (x, z) def = j ∈ [1, n i ]|C ⊤ j x + µ k z k j + µ k (α -1)z j -u j > 0
refers to the primal-dual active set at (x, z) for the current k th sub-problem, and J k (x, z) c stands for its complementary set (i.e., the set of inactive constraints). The linear system ( 14) is symmetric and non-singular. Thus it can be efficiently solved using a Cholesky factorization. Furthermore, one can notice it does not explicitly contain any square matrix multiplications. We will see it is not the case when working with Φ k ρ,µ k . At the l th iteration, a semi-smooth Newton method applied to Φ k ρ,µ k , and initialized at

x(0) = x k involves finding dx satisfying ∇ 2 Φ k ρ,µ k (x (l) )dx + Φ k ρ,µ k (x (l) ) = 0.
One ends-up solving [START_REF]QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs[END_REF]Equation 3.4]:

H + ρI C ⊤ I k (x (l) ) C I k (x (l) ) -µ k I dx dz = -∇ x Φ k ρ,µ k (x (l) ) 0 ,
where

I k (x (l) ) def = {i ∈ [1, n i ] | C i x(l) -u i + z k i µ k > 0}
refers to the active set of the current proximal sub-problem, and C I k (x (l) ) corresponds to a reduced version of the matrix C containing the active rows indexed by I k (x (l) ). The gradient on the right-hand side of this equation equals:

∇ x Φ k ρ,µ k (x l ) =H x(l) + g + ρ(x (l) -x k ) + C ⊤ z k + 1 µ k (C x(l) -u) + .
For either low values of µ k or ill-conditioning of C I k (x (l) ) , the evaluation of the gradient is numerically challenged due to the value of

1 µ k C ⊤ I k (x (l) ) C I k (x (l)
) (see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 17.4] which gives examples illustrating ill-conditioning resulting from such structures).

In the next section, we present other state-of-the-art methods for solving QPs alongside their theoretical features that can be used in the context of embedded applications.

IV. RELATED WORK

Generic-constrained convex QPs are commonly solved with iterative (or indirect) methods. These methods are traditionally divided into two main families: (i) active-set methods, and (ii) penalization methods. When the problem contains no inequality constraint (i.e., unconstrained or only subject to equality constraints), direct methods can be used, as solving the corresponding QP can be done by directly solving the corresponding linear KKT system [START_REF] Nocedal | Numerical optimization[END_REF]. Iterative methods of both families typically aim at solving cascades of structurally simpler intermediary optimization problems, whose solutions eventually go to that of the original problem.

A. Active-set methods

This family of methods, first developed in the 50s, aims at determining the set of active constraints at an optimal point of (QP) (see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.5] for an introduction). Once this set of active constraints is determined, the constrained QP can be solved directly by considering an equivalent QP subject to equality constraints only, whose solution matches that of the original QP. On the positive side, warm-start and hot-start strategies can easily be incorporated within activeset methods [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF]. This is usually key for solving cascades of similar QPs, which is common in applications such as sequential quadratic programming (SQP) and model predictive control (MPC). They also usually perform well on small-sized QPs. Yet, active-set methods have limited scalability, struggle to exploit sparsity, and are hard to early-stop at inexact solutions satisfying (1) at accuracy ϵ abs (they tend to stop only when a global solution is found, i.e., with ϵ abs = 0). Moreover, they may require strong assumptions, such as the QP being strongly convex or satisfying "linear independence constraint qualification" (LICQ) (see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Definition 12.4]). These assumptions reduce their applicability and may cause robustness issues, such as "active-set cycling" [28, sections 13.5 and 16.5].

B. Penalization methods

This family of methods transforms the original constrained problem into a sequence of problems with either no constraints or very simple ones (such as sign constraints on some variables). The new problems typically have objectives that consist of two terms: (i) the objective of the original problem and (ii) a term for penalizing points not being feasible for the original problem. There are two very common types of penalization methods used in practice.

Interior-point methods. This family of penalization methods forces through a barrier function the sequence of intermediary optimization problems to have strictly feasible solutions with respect to the domain of the original problem; see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.6] for an introduction. Primal-dual interior-point methods [START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF], [START_REF] Wright | Primal-dual interior-point methods[END_REF] became popular in the 90s due to their good practical performances across a wide range of problems. On the positive side, unlike active-set methods, they benefit from early termination and can easily exploit the sparsity structure of the problem. Moreover, by construction, they control the primal-dual gap. Yet, because of the homotopy-based structure of this family of methods, one of their main drawbacks is the difficulty of using them with warm-starting procedures when solving a sequence of related QPs (e.g., within SQP solvers or MPC problems). Augmented Lagrangian-based methods. This family of methods is primarily based on the idea of Lagrangian relaxations with an additional quadratic penalization term (possibly piecewise) for encouraging the feasibility of the iterate (see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Section 17.3]). This kind of technique emerged in the 70s through the works of Hestenes and Powell [START_REF] Hestenes | Multiplier and gradient methods[END_REF], [START_REF] Powell | A method for nonlinear constraints in minimization problems[END_REF] and then later with those of Rockafellar [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF], which tightly emphasized their link with the so-called "proximal-point method". Indeed, augmented Lagrangian-type methods naturally arise by applying proximal-point methods on either the dual or saddlepoint formulations to (QP), thereby offering the advantage of converging under relatively weak assumptions. Augmented Lagrangian-type methods also have the advantage of exploiting the sparsity structure of the problem and can be easily warmstarted via warm-start and hot-start procedures. They also benefit from early termination. Moreover, as shown in Remark 1 previously, augmented Lagrangian algorithms based on the method of multipliers inherit from the advantageous property of converging towards the closest feasible QP in ℓ 2 sense, if it happens to be primal infeasible [START_REF] Chiche | How the augmented Lagrangian algorithm can deal with an infeasible convex quadratic optimization problem[END_REF]. However, AL-based methods might exhibit slow convergence behaviors similar to those of first-order optimization methods. For instance, when they are based on splitting methods such as the alternating direction method of multipliers (ADMM) [START_REF] Parikh | Proximal algorithms[END_REF]Section 4.4] they are relatively slow toward high accuracy levels. Furthermore, they might not necessarily have control over the primal-dual gap.

V. SOFTWARE IMPLEMENTATION

In this section, we list a range of common QP solvers, along with their underlying optimization strategy, then detail the software implementation of PROXQP within the PROXSUITE library.

A. Common QP solvers

We provide below a summary of different start-of-the-art solvers implementing iterative optimization methods previously described in Section IV, and summarized in Table I, with a short description of their practical features for embedded optimization.

Active-set solvers. Popular active set-based convex QP solvers include the open-source QPOASES [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF], QUADPROG [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF], DAQP [START_REF] Arnström | A dual active-set solver for embedded quadratic programming using recursive LDL T updates[END_REF], and the QPA module in the open source software GALAHAD [START_REF] Gould | GALAHAD user documentation[END_REF]. Interior-point solvers. Standard solvers using an interiorpoint strategy are commercial solvers GUROBI [START_REF] Optimization | GUROBI optimizer reference manual[END_REF] and MOSEK [START_REF] Mosek | MOSEK optserver documentation[END_REF], closed-source BPMPD [START_REF] Mészáros | The bpmpd interior point solver for convex quadratic problems[END_REF], and open-source solvers OOQP [START_REF] Gertz | Object-oriented software for quadratic programming[END_REF], HPIPM [START_REF] Frison | Hpipm: a high-performance quadratic programming framework for model predictive control[END_REF], ECOS [START_REF] Domahidi | ECOS: an SOCP solver for embedded systems[END_REF], CVXOPT [START_REF] Andersen | Cvxopt: A python package for convex optimization[END_REF] and QPSWIFT [START_REF] Pandala | qpSWIFT: A Real-Time Sparse Quadratic Program Solver for Robotic Applications[END_REF]. Augmented Lagrangian-based solvers. Common solvers based on augmented Lagrangians for solving QPs include OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF] and SCS [START_REF] O'donoghue | Conic optimization via operator splitting and homogeneous self-dual embedding[END_REF] (via an alternating direction method of multipliers) and QPALM [START_REF]QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs[END_REF].

B. The PROXSUITE library

PROXQP is implemented in C++ within the PROXSUITE library and builds extensively on the Eigen library [54] for efficient linear algebra. The current implementation is tailored to both dense and sparse matrix operations and leverages recent CPU architectures equipped with advanced vectorization instructions. PROXSUITE is publicly available at https://github.com/Simple-Robotics/proxsuite under a permissive open source BSD-2-Clause license. It comes with out-ofthe-box interfaces for the C++, Julia, and Python languages inspired from OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]. The documentation of PROXSUITE's API is available online at https://simple-robotics.github.io/ proxsuite/. In the following, we describe some key features of PROXSUITE, list our default parameters in Table II as well as a code snippet in Listing 1. Cholesky factorization. In PROXSUITE, we have developed dedicated dense and sparse LDLT Cholesky factorizations [START_REF] Golub | Matrix computations[END_REF] to explicitly account for the specific calculations. In particular, this Cholesky factorization implements advanced rank update routines to efficiently account for the change of active sets when solving [START_REF] Mayne | Robust model predictive control of constrained linear systems with bounded disturbances[END_REF] while lowering the overall memory footprint. Furthermore, when the number of constraints is too large compared to the number of (primal) variables, we add an option to factorize, using the dense backend, the matrix

H + ρI + 1 µ C ⊤ J k (x (l) ,ẑ (l) ) C J k (x (l) ,ẑ (l) ) , (15)
instead of the one involved in ( 14) (the latter choice is called PRIMALDUALLDLT, whereas the first PRIMALLDLT).

A heuristic (comparing "typical" algorithmic complexity for factorizing, updating, and solving (14) using one of the mentioned matrices) is also available for automatically choosing which matrix to factorize. Moreover, a similar heuristic is available for the sparse backend when the system is large-scale. In such case, a sparse matrix-free routine then solves [START_REF] Villa | Model predictive control of biped walking with bounded uncertainties[END_REF] in the same spirit as [START_REF]QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs[END_REF]Section 4.1.3]. The sparse and dense Cholesky factorizations are also freely available within the PROXSUITE library and will hopefully be integrated into the Eigen library soon. Warm-start. Since PROXQP is built on an augmented Lagrangian-based method, the algorithm is readily warmstartable. For example, it is a useful feature for solving a cascade of QPs that change slightly (e.g., for Model Predictive

[40] Dense ✓ ✓ ✗ ✓ ✗ ACTIVE-SET QUADPROG [44] Dense ✗ 1 ✗ 1 ✗ ✓ ✗ DAQP [45] Dense ✗ 1 ✗ 1 ✗ ✓ ✗ GUROBI [22] Sparse ✗ ✗ ✓ ✓ ✗ INTERIOR- MOSEK [21] Sparse ✗ ✗ ✓ ✓ ✗ POINT CVXOPT [51] Dense ✗ ✗ ✓ ✓ ✗ ECOS [50] Sparse ✗ ✗ ✓ ✓ ✗ QPSWIFT [52] Sparse ✗ ✗ ✓ ✓ ✗ OSQP [20] Sparse ✓ ✓ ✓ ✗ ✗ 2 AUGMENTED SCS [53] Sparse ✓ ✓ ✓ ✓ ✗ 2 LAGRANGIAN PROXQP [16] Sparse & Dense ✓ ✓ ✓ ✓ ✓
1 In their original implementation. 2 Not established nor implemented. penalization for inequality constraints Control, Inverse Kinematics, Inverse Dynamics, SQPs, etc.). Default initial guesses are also available as options. Among others, by default and motivated by the fact that equality constraints are always active at an optimal solution, we use the following initialization for primal and dual variables:

x 0 z 0 = H + ρI C ⊤ C -µ 0 I -1 -g u ,
provided C is an equality constraint matrix and u and equality constraint vector. Hot-start. PROXQP has a dedicated feature, similar to the OSQP, SCS, and QPOASES solvers [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF], [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF], [START_REF] O'donoghue | Conic optimization via operator splitting and homogeneous self-dual embedding[END_REF], where we can hot-start the solver with a factorization from a previous related problem. More precisely, if QP updates concern only linear terms (e.g., u or g), then the updates are realized allocation-free in PROXQP. Furthermore, PROXQP has settings that combine warm-starts with a reuse of the whole previous data workspace (including the Cholesky factorization) called WARM START WITH PREVIOUS RESULT. This feature has been specially designed for solving control or SQP problems with changes concerning only linear parts.

Preconditioning. PROXQP contains several preconditioning strategies, enhancing the overall numerical stability and convergence of the optimization process. The default preconditioner used in our current implementation is often referred to as the Ruiz equilibration [START_REF] Ruiz | A scaling algorithm to equilibrate both rows and columns norms in matrices[END_REF]; see, e.g., [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]Algorithm 2], also used in OSQP and QPALM.

Solving closest feasible problems. PROXQP can solve the closest feasible QP problem if it is detected to be primal infeasible (following the same procedure described as in OSQP [20, Section 3.4]). If the option is activated, then the following stopping condition will be tracked:

   ∥Hx ⋆ + g + C ⊤ z ⋆ ∥ ∞ ⩽ ϵ abs ∥C ⊤ [Cx ⋆ -u] + ∥ ∞ ⩽ ∥C ⊤ 1∥ ∞ ϵ abs |(x ⋆ ) ⊤ Hx ⋆ + g ⊤ x ⋆ + (u + s ⋆ ) ⊤ [z ⋆ ] + | ⩽

VI. EXPERIMENTAL RESULTS

In this section, we evaluate PROXQP against other stateof-the-art approaches on classic control problems and representative robotic tasks. These practical problems illustrate how PROXQP takes advantage of structural properties (sparsity structure, warm-starts, hot-starts, early stopping). We also show that PROXQP's ability to solve primal infeasible problems can be leveraged in closed-loop MPC scenarios accounting for perturbations and measurement errors that yield infeasible problems. Finally, we conclude benchmarking PROXQP against alternative solvers of the state of the art on more generic types of QPs [START_REF] Maros | A repository of convex quadratic programming problems[END_REF], [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF].

A. Benchmark setup and scenarios

Our benchmarks focus on relatively small to mediumsized QPs, typically used for embedded applications such as in robotics (less than a thousand variables and constraints), with sparse or dense problem structures.

In the first set of experiments (see Section VI-B), we consider sparse convex MPC scenarios: the benchmark proposed by [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]Section 8.4] and the one proposed by [58, Section 5.A]. We also explain with these synthetic examples how primal infeasibility solving can be used in the context of closed-loop MPC : (i) first, we initialize one instance with a primal infeasible point with respect to the dynamics bounds x min and x max ; (ii) in the latter, we add an infeasible random perturbation to the dynamics at a fourth of the simulation.

In the second set of experiments, we focus on dense problems (see Section VI-C): inverse kinematics tasks, a closed-loop dense MPC tasks for controlling the center of mass of a humanoid robot, and a real-world closed-loop MPC scenario with the Upkie wheeled-biped robot [START_REF]Upkie wheeled biped robots: code and build instructions[END_REF]. Through an extensive benchmark with the humanoid closed-loop MPC task, we also highlight that PROXQP is more robust than other state-ofthe-art solvers against random perturbations that render the approximate solution infeasible.

Finally, we evaluate PROXQP on more generic and standard QPs of the literature from the optimization community.

For all experiments, PROXQP is compared against other state-of-the-art approaches: active-set methods (QPOASES, QUADPROG), interior-point methods (MOSEK, QPSWIFT) and augmented Lagrangian-based methods (SCS, OSQP). We exploit all available features (see Table I) of each solver considering the problem structure. More precisely, OSQP, QPOASES, and SCS use warm-starting and hot-starting when possible. Box constraints on the state or control variables are specified when the solver API allows to do so. Otherwise, they are listed as parts of the inequality constraints. Moreover, interior-point and augmented Lagrangian-based methods exploit their underlying early-stopping strategy if it suits the task. The experiments were conducted on a standard laptop CPU, an Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz. The benchmark is open-source and available at https://github.com/Bambade/proxqp benchmark with an easy-to-use interface strongly inspired by the one proposed by OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF].

B. Sparse problems

In this subsection, we consider two sets of convex MPC problems: an existing convex-MPC benchmark from [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]Section 8.4], and a sparse chain of oscillating masses from [START_REF] Wang | Fast model predictive control using online optimization[END_REF]Section 5A]. In both cases, we run systems over 100 simulation steps. Their initial states x init are uniformly chosen between reasonable bounds. We thus measure the average of the full-time needed for solving all 100 simulation steps. For each separate randomization seed, the full simulation steps' is also averaged over 10 consecutive runs in order to mitigate the impact of loading costs. Furthermore, we consider two different scenarios: (i) one when early-stopping is relevant (ϵ abs = 10 -3 ); and another (ii) where higher accuracy is required (ϵ abs = 10 -6 ). In this test set, only the vectors terms in (QP) change in this test case, the practitioner can use warm-start and hot-start strategies if available. These features are exploited by OSQP, QPOASES, SCS and PROXQP. Furthermore, inequality constraints in these problems only consist of box constraints, a specificity that SCS leverages with a custom feature. Finally, since the problems have a sparse structure, PROXQP uses its sparse back-end.

1) Convex MPC test set Experimental setup. In this set of experiments, we consider the MPC benchmark proposed in [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]Section 8.4]. This benchmark consists of LQR problems with additional state and input constraints, as described in [20, Appendix A3], with a horizon T = 10 starting from

x 0 = x init ∼ U(-0.5x, 0.5x) with x ∼ U(1, 2): min xt∈R nx ,ut∈R nu x ⊤ T Q T x T + T -1 t=0 x ⊤ t Qx t + u ⊤ t Ru t , s.t. x t+1 = Ax t + Bu t , -x ⩽ x t ⩽ x, -ū ⩽ u t ⩽ ū, (16)
with

Q, Q T ∈ S + (R nx ), R ∈ S ++ (R nu ), and A ∈ R nx×nx , B ∈ R nx×nu .
More precisely, we solve this optimal control problem in a receding horizon fashion, computing an optimal input sequence u 0 , . . . , u T -1 , applying the first input u 0 to the system and propagating the state to the next time step. The whole procedure is repeated 100 times. We choose the state variable dimension n x and control input dimension n u for the context of small systems embedded systems, concretely n x = 6, 12, 16, 20 and n u = n x /4. This benchmark generates sparse problems (i.e., the KKT matrix of the QP has about 5% of nonzeros elements) and QPs whose dimensions (number of primal variables and constraints) range from 76 to 490. Runtime results. The first and second-row blocks of Table III show the results of this test case. We can see that in both scenarios PROXQP is the fasted method with a speed-up ranging from 1.2 to 1.8 in the early-stopping scenario (with respect to OSQP, the second fastest approach); and from 1.7 to 2.9 in the high accuracy scenario. We can also observe that sparse solvers (PROXQP, OSQP, SCS, QPSWIFT, MOSEK) obviously provide better timings than dense solver (QPOASES). Furthermore, one can notice that methods based on ADMM (i.e., OSQP and SCS) are slower to converge to high-accuracy solutions. 1 The solution does not satisfy the desired absolute accuracy.

2 QUADPROG cannot solve QPs that are not strictly convex.

Solving MPC from primal infeasible points. In order to solve an MPC problem starting from x init ∼ x + U(0, 0.01), we consider the following strategy: the solution found at the first simulation step provides a solution trajectory xt , ût (over a horizon of T time steps for control inputs ût and T + 1 for the state xt ) for the closest feasible QP of ( 16) (see orange line in Figure 1). PROXQP also provides optimal shifts ŝu t , ŝx t for the inequality constraints on x and u; and ŝe t for equality constraints. We make the forward updates as follows:

• u 0 = û0 -ŝu 0 : we bring back the closest feasible input solution to the solution space of the initial QP (so u 0 satisfies the input constraints since û0 ⩽ ū + ŝu 0 by construction); • x 1 = Ax 0 +Bu 0 -ŝe 0 : we update the trajectory from a feasible control input u 0 and optimal dynamic x0 and shift ŝe 0 . It lies on the green line drawn in Figure 1. Compared to the closest feasible optimal solution x1 = Ax 0 + B û0 -ŝe 0 , we make a drift

x 1 -x1 = Ax 0 + Bu 0 -ŝe 0 -[Ax 0 + B û0 -ŝe 0 ], = B(u 0 -û0 ) = -Bŝ u 0 ,
which tends towards zero, as soon as the optimal shift over the control input tends towards zero. • We repeat the closed-loop process starting from x 1 to derive x1 , û1 and then

x 2 = Ax 1 + B(û 1 -ŝu 1 ) -ŝe 1 .
We illustrate the application of this update rule with an instance of the convex MPC scenario for n x = 6 and n u = 1 in Figure 2. Figure 2b shows the control trajectory u t of the system. As expected, it lies in the control bounds -ū and ū. Figure 2a shows the 6 components of the dynamics x t . All these components are initially infeasible for a fixed number of simulation steps since they are above x. Since the initial infeasible perturbation is not too large (i.e., the optimal shifts are sufficiently small), the dynamics eventually converge to the feasible steady-state position.

2) Sparse chain of oscillating masses Experiment setup. We consider the following MPC problem involving a chain of six masses interconnected with spring- x 1 ~time step state u Fig. 1: A convex MPC problem is initialized with an infeasible state x 0 = x init > x max . Grey: Yet, x LQR 0 = x 0 enables reaching a horizon T a feasible state x goal with an unconstrained LQR scheme (see (??)). Orange: It guarantees that a closest feasible trajectory xt exists with optimal shifts ŝx t ,ŝ u t ,ŝ e t . Green: We choose the update x 1 def = x1 -Bŝ u 0 , which is guaranteed not moving away x 1 from x goal for sufficiently small optimal shifts ŝx 1 ,ŝ u 0 ,ŝ e 1 . More precisely, ∥x 1 - Appendix-B) provides more technical details about the properties of this update rule. dampers [58, Section 5A] III show the runtime results. We can see that PROXQP is the fastest method with a speed-up ranging from 1.3 (early stopping scenario) to 1.5 (high accuracy scenario).

x goal ∥ 2 ⩽ ∥x LQR 1 -x goal ∥ 2 and ∥x 1 -x max ∥ 2 ⩽ ∥x LQR 1 -x goal ∥ 2 . (
min x t ∈R nx ,u t ∈R nu 1 2 x ⊤ T Q T x T + 1 2 T -1 t=0 x ⊤ t Qxt + u ⊤ t Rut, s.t. x t+1 = Axt + But, x 0 = x init ∼ U (-0.5, 0.5), - 4
Solving MPC with infeasible perturbations. To illustrate another possible scenario where primal infeasibility solving might be useful, we consider the following setup: at bout one fourth of the simulation (i.e., the 25 th simulation step), a random perturbation w t ∼ 4.8 + U(0, 0.01) disturbs the usual and n u = 1 starting from an infeasible initial state x 0 using ProxQP. At each new simulation step t, the new control input equals u t = ût -ŝu t , so it is guaranteed to be feasible, while the dynamics follow x t = xt -Bŝ u t . Since the initial infeasible perturbation is not too large, the dynamics converge after a fixed number of steps to a feasible regime. dynamic update

x t+1 = Ax t + Bu t + w t ,
which is thus guaranteed at the next LQR simulation step to create an infeasible problem since x t is already close to the steady state for such a system after 25 steps. We apply the same update rule as exposed in the previous MPC problem to tackle this infeasible perturbation since the problem structure is similar. Figure 3a shows the 12 components of the chain system. It becomes infeasible after the 25 th closed-loop update. Nevertheless, we can see that it eventually recovers to the steady state solution after a finite number of time steps. Figure 3b shows the 3 control input components. We can see that, as expected, they all satisfy the constraint bounds thanks to the update rule u t = ût -ŝu t .

C. Dense problems 1) Inverse kinematics

In this set of experiments, we consider different inverse kinematic (IK) tasks proposed for motion control of articulated robots using Pinocchio [START_REF] Carpentier | The pinocchio c++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives[END_REF] and the open-source Pink library https://github.com/tasts-robots/pink. The resulting motion controller is run for 5 s (i.e., solving 1, 000 consecutive QPs Fig. 3: Convex MPC problem for controlling a chain of masses with an infeasible dynamic perturbation at the fourth of the simulation. The control inputs and the dynamics are updated following paragraph VI-B1c. After a fixed number of time steps, the system manages to recover its asymptotic steady state while satisfying during the whole process feasible control bounds.

with a time-step of 5 ms). Each such run is repeated 50 times, and we measure the runtime per timestep average for achieving each IK task. Considering the structure of the problem, OSQP, SCS, QPOASES, and PROXQP benefit from warm-starting, and PROXQP uses its dense-backend. Furthermore, from a practitioner's point of view, we consider the minimal accuracy needed ϵ abs for achieving these tasks (i.e., solving each IK task while maintaining sufficient smoothness in the resulting motions). We thus imposed ϵ abs = 10 -5 for each task, except when using STRETCH which required a lower ϵ abs = 10 -7 in order to satisfy to be numerically stable while enforcing joint limits.

Table IV reports runtimes in its second-row block. It shows that PROXQP is the fastest method with a speed-up ranging from 1.2 to 1.8 with respect to OSQP or QUADPROG, the next two fastest approaches. Apart from OSQP, we note how dense solvers are faster than sparse ones, which can be explained by the fact that IK problems depict a relatively high ratio of nonzero elements with respect to the problem size (the ratios of sparsity of the Hessian cost matrix H ranges from 38% to 100%, and between 4% and 16% for the constraint matrix).

2) Chain of oscillating masses

In this set of experiments, we consider a dense version of the chain of oscillating masses MPC problem proposed by [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF]Section 6.3] as an SQP problem. The first test problem CHAIN80 aims at regulating a chain of nine masses connected by springs into a certain steady state. One end of the 1 The solver throws a factorization error (of non-convexity). 2 The solution does not match the desired accuracy 3 The solver does not manage to satisfy configuration limits. 4 The solver does not manage to handle upper bound constraints and outputs infeasibility errors. 5 Low accuracy iterates provokes with SCS warm starts more difficult QPs to solve in a closed-loop strategy. QPs with 240 variables and 709 constraints. These two SQPs are run 10 times, and we measure the total runtime necessary for solving these tasks to an accuracy ϵ abs = 10 -6 . We can see the runtime results in the first-row block of Table IV. It shows that PROXQP is the fastest method with a speed-up ranging from 1.8 to 11.6 with respect to QPOASES, the second fastest approach. We can notice again that dense solvers are faster than the other sparse solvers.

3) Closed-loop MPC

In this last set of experiments, we consider dense convex MPC problems for controlling a walking humanoid robot in simulations, and the real robot Upkie [START_REF]Upkie wheeled biped robots: code and build instructions[END_REF] in a balancing scenario. We evaluate again the robustness of the solvers to perturbations.

Humanoid locomotion via cart-table model. The task consists in controlling a reduced model of a humanoid robot, consisting of its center-of-mass c and zero-tilting moment point (ZMP) z, using the center-of-mass jerk u = ... c as control input [START_REF] Kajita | The 3d linear inverted pendulum mode: a simple modeling for a biped walking pattern generation[END_REF], [START_REF] Caron | Stair climbing stabilization of the hrp-4 humanoid robot using whole-body admittance control[END_REF]. Figure 4 depicts a trajectory of this system. The goal is to keep the ZMP inside the support polygon of the feet while achieving a prescribed walking velocity. The resulting MPC problem can be cast as a quadratic program [START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF] of the form:

min x t ∈R nx u t ∈R nu wT ∥xT -xgoal∥ 2 2 + T -1 t=0 wx∥xt -xgoal∥ 2 2 + wu∥ut∥ 2 2 , s.t. xt+1 = Axt + But, x0 = xinit, Cxt + Dut ⩽ et, (17)
where n x = 3 (for center-of-mass position c, velocity ċ and acceleration c), n u = 1 and T = 16. As in [START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF], horizontal coordinates can be decoupled, so that we can focus on singledimensional coordinates c, z ∈ R. For the sake of this example we choose lateral coordinates and a target velocity of zero corresponding to walking in place.

The problem is solved with respect to the stacked vector U = [u 0 . . . u T -1 ] of control inputs, deriving states x t linearly from U by recursive expansion of the dynamics x t+1 = Ax t + Bu t [START_REF] Audren | Model preview control in multi-contact motion-application to a humanoid robot[END_REF]. This leads to a series of dense QPs with 16 variables and inequality constraints. The experiment is run for 100 simulation steps with a duration of 0.1s each, and the results are averaged over 100 trials. We measure the average total runtime for solving the task. Since the system is time-invariant, OSQP, SCS, QPOASES, and PROXQP benefit from warm-starting and hot-starting. In our experiments, setting ϵ abs = 10 -2 is sufficient to solve the task for all solvers based on IP and AL-based methods. Yet, they do not benefit the same way from this early-stopping feature. For this reason, we execute the runtime benchmark for ϵ abs = 10 -2 , ϵ abs = 10 -4 and ϵ abs = 10 -6 .

The runtime results are displayed in the last row block of Table IV. PROXQP is slightly faster than OSQP or QUADPROG, the best performing state-of-the-art solvers on this test set, where we note again how, as expected, dense solvers are generally faster than sparse solvers. Finally, we can see that SCS benefits less from early stopping here. Timings for SCS improve when the iterates are more accurate. Robustness to perturbations. We apply random (Gaussian) noise at each simulation step to the center-of-mass acceleration with an amplitude σ modeling state observation inaccuracies or unmodeled dynamics. When the ZMP is close to the edge of its support polygon, a small perturbation in center-of-mass acceleration can cause it to cross the edge, which would break the foot contact. For this reason, it is common practice to add safety margins to ZMP support areas in real systems, sometimes reducing them to a square well within the full area [START_REF] Scianca | Mpc for humanoid gait generation: Stability and feasibility[END_REF]. We include such margins in our QP formulation, as depicted by the light-dotted lines in Figure 4. In order to quantify the robustness of the different solvers to various noise amplitudes, we measure the percentage of the full trajectory (100 simulation steps) each solver manages to solve and average over 100 seeds. This closed-loop task is interrupted if the primal feasibility violation with respect to the original support polygon (dark dotted lines) exceeds the accuracy of ϵ abs = 10 -2 , or if the solver does not return any solution. For all solvers, we set the primal infeasibility tolerance to the minimum possible value in order to solve as many problems as possible. In Table V Fig. 5: Closed-loop model predictive control with PROXQP on the Upkie wheeled biped. The top plot shows the time taken by the hot-started solver to solve MPC problems at each control cycle on the robot's embedded computer. The bottom plot shows the angle between the robot's torso and the vertical at the corresponding times, illustrating how it balances upright. Upkie wheeled biped. The controller runs entirely on the robot's embedded computer, a Raspberry Pi 4 Model B equipped with a quad-core ARM Cortex-A72 64-bit SoC @ 1.8 GHz. Note that this CPU has fewer cores and a lower frequency than the standard laptop CPU used in other test cases. We formulate the MPC problem as a quadratic program using the same strategy as in [START_REF] Bambade | QPLayer: efficient differentiation of convex quadratic optimization[END_REF], for a discretization of the linearized cart-pole model (for wheels) rather than cart-table (for flat feet). States observed from IMU Kalman filtering and wheel odometry are fed directly as an initial state x init to the MPC problem. We set T = 50 timesteps with a time discretization of 20 ms, resulting in a receding horizon duration of 1 s. Task weights are chosen as w T = 1, w x = 10 -3 and w u = 10 -3 . The resulting QPs have dimension n = 50 with m = 100 inequality constraints. They are solved by PROXQP with hot-starting in 0.8±0.02 ms, fitting well within the budget of a 200 Hz control loop. Figure 5 reports the resulting computation times alongside observed states for a sample run of the balancing controller on the robot. The source code for this experiment is available in the mpc_balancer directory of the robot's repository [START_REF]Upkie wheeled biped robots: code and build instructions[END_REF].

D. Generic QPs 1) Random QPs

We generated random QPs with 15% of sparsity for both matrices H and C in the spirit of the benchmark API from OSQP [START_REF] Bambade | PROX-QP: Yet another Quadratic Programming Solver for Robotics and beyond[END_REF], [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF], [START_REF]QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs[END_REF]. The dimension d of those problems is in the range from 10 to 1000. For each dimension, we 6: Performance profiles on small to medium-sized Maros-Mészáros problems, using a Core i5 -5300U -5th Generation @ 2,3 GHz processor. Target accuracy ϵ abs = 10 -6 and time limit set to 1000 seconds. The higher the better. generated 1000 problem instances with different random seeds and averaged the runtime of each algorithm on 10 consecutive runs to limit the impact of loading costs. We then report the failure rate and the shifted geometric means [16, Section VI-Db] in Table VI. PROXQP is about 2.6 times faster across all the problems than OSQP, the second fastest approach. Furthermore, PROXQP, OSQP, and QPOASES exhibit no failure rate.

2

) Maros-Mészáros problems

The Maros-Mészáros test set [START_REF] Maros | A repository of convex quadratic programming problems[END_REF] is composed of 138 "hard" QPs. Most of them are sparse and ill-conditioned problems, and they contain up to 90597 variables and 180895 constraints. About 83% of the problems have a sparsity level lower than 10% for the Hessians H, as well as for the constraint matrices. We restrict the benchmark to problems whose dimensions (constraints and variables) are below or equal to 10 3 . This results in a subset of 62 problems, about 45% of the Maros-Mészáros test set, for which two-third of the problems have a sparsity level no larger than 20%.

We report performance profiles [START_REF] Bambade | PROX-QP: Yet another Quadratic Programming Solver for Robotics and beyond[END_REF], [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF], [START_REF]QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs[END_REF] on Figure 6 comparing PROXQP with its dense backend against the other solvers of the testbed. Performance profiles correspond to the fraction of problems solved as a function of a certain runtime (measured in terms of a multiple of the runtime of the fastest solver for that problem). Figure 6 depicts that even-though the problems are sparse, PROXQP always depicts the best performance profile. Hence it is the fastest approach for solving the problems of this testbed.

VII. CONCLUSION

In this work, we have introduced a new algorithm for solving generic QPs, motivated by robotic and control applications.

Among others, we combine the bound constraint Lagrangian (BCL) globalization strategy [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF] for automatically scheduling key parameters of a primal-dual proximal augmented Lagrangian algorithm. The intermediary proximal subproblems are solved via a semi-smooth Newton method. The convergence guarantees come together with an efficient software implementation of PROXQP in the open-source PROXSUITE library, which contains several practical features for embedded optimization. The performance of PROXQP is evaluated on different standard sparse and dense robotic and control experiments, including a real-world closed-loop controller application. In particular, we have highlighted how solving the closest feasible QPs can be efficiently leveraged in the context of closed-loop convex MPC. Finally, through various benchmarks of the optimization literature [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF], [START_REF] Maros | A repository of convex quadratic programming problems[END_REF], [START_REF] Wang | Fast model predictive control using online optimization[END_REF], we have shown that PROXQP performs at the level of state-of-the-art solvers on a large set of generic QP problems.

As future works, we plan to further study and extend the primal-dual proximal augmented Lagrangian approach developed in this work to more general contexts beyond quadratic costs and linear constraints, as initiated through recent works [START_REF] Jallet | Constrained differential dynamic programming: A primal-dual augmented lagrangian approach[END_REF], [START_REF] Jallet | Proxnlp: a primal-dual augmented lagrangian solver for nonlinear programming in robotics and beyond[END_REF]. Another appealing research direction concerns the extension of the approach for proper differentiation of QP problems in the context of differentiable optimization [START_REF] Agrawal | Differentiable convex optimization layers[END_REF], notably in the context of infeasible QPs [START_REF] Bambade | QPLayer: efficient differentiation of convex quadratic optimization[END_REF]. From a software perspective, we plan to provide long-term support to help disseminate the approach for research, industrial, and commercial contexts.

Algorithm 1 :

 1 PROXQP (practical implementation)

⩽ xt ⩽ 4 ,

 4 -0.5 ⩽ ut ⩽ 0.5, with R = I and Q = I. The problem dimensions are n x = 12, n u = 3 and the horizon is T = 30 [58, Section 5A]. It generates QPs with 462 variables and 864 constraints. Runtime results. The last two-row blocks of Table

6

 6 components of the dynamics x t over 100 simulating steps.

  Control input u t over 100 simulating steps.

Fig. 2 :

 2 Fig.2: Closed-loop convex MPC simulation steps with n x = 6 and n u = 1 starting from an infeasible initial state x 0 using ProxQP. At each new simulation step t, the new control input equals u t = ût -ŝu t , so it is guaranteed to be feasible, while the dynamics follow x t = xt -Bŝ u t . Since the initial infeasible perturbation is not too large, the dynamics converge after a fixed number of steps to a feasible regime.

  State dynamic components over 100 simulation steps.

  Control input components over 100 simulation steps.

  u max u min (b) 3 components of the inputs u t over 100 simulation steps.

Fig. 4 :

 4 Fig.4: Controlling the lateral center-of-mass trajectory (blue) to maintain the ZMP (red) within the real support polygon (dotted dark green and blue). Light dotted lines are more conservative ZMP bounds used during solving QP. Random perturbation of the current center-of-mass acceleration can lead to a ZMP outside the support polygon and a primal infeasible QP.

TABLE I :

 I Comparison of different optimization methods and software libraries for solving QPs.

	Method	Solver	Backend	Warm-	Hot-starting	Early	Primal-dual	Infeasibility
				starting		termination	gap	solving
		QPOASES

TABLE II :

 II Default hyperparameters for PROXQP.

	Parameter	Value	Description
	α	0.95	primal-dual interpolation
	ϵ 0 bcl = ϵ 0 ρ	1 10 -6	accuracy primal proximal penalization
	µ min	10 -9	minimal penalization
	µ f	0.1	penalization decay factor
	α bcl , β bcl kmax	0.09, 0.9 10 4	BCL hyperparameters safeguard
	µ 0 e µ 0 i	10 -3 10 -1	penalization for equality constraints

TABLE III :

 III Sparse convex MPC benchmark (total runtimes in ms for solving 100 simulation steps)

	Problem type	Dimensions	PROXQP	QUADPROG	OSQP	QPOASES	SCS	QPSWIFT	MOSEK
		n=76,m=142	5.7±0.2	✗ 2	6.7±0.7	37.6±23.9	11.9±2.2	13.5±0.5	82.1±3.5 1
	CONVEX MPC	n=162,m=294	16.4±1.7	✗ 2	24.1±6.0	393.9±109.9	51.0±15.7	47.1±4.6	211.5±13.3
	ϵ abs = 10 -3	n=216,m=392	27.6±4.8	✗ 2	48.2±29.4	1434.2±480	65.6±30.9	80.0±3.6	311.3±15.4 1
		n=270,m=490	43.2±8.3	✗ 2	77.5±65.6	2345.5±608	104.7±66.0	131.0±7.4	454.2±18.6
		n=76,m=142	8.5±1.4	✗ 2	14.1±8.4	27.1±15.4	25.5±14.6	16.5±0.7	84.0±3.8 1
	CONVEX MPC	n=162,m=294	21.1±2.3	✗ 2	94.0±101.9	452.9±98.3	101.8±47.9	57.8±5.3	216.3±14.1 1
	ϵ abs = 10 -6	n=216,m=392	36.1±4.3	✗ 2	150.1±136.0	1626.6±650.1	162.0±95.0	99.5±4.9	314.1±17.4 1
		n=270,m=490	57.0±9.0	✗ 2	346.7±362.0	2670.3±842.4	311.1±122.5	164.1±9.9	459.6±21.0 1
	CHAIN OF MASS	n=462,m=834	45.2±0.8	(28.3±0.6)E3	59.1±1.7	(80.2±1.6)E3	161.1±11.1	215.3±1.3	566.8±14.1
	ϵ abs = 10 -3
	CHAIN OF MASS	n=462,m=834	67.8±11.1	(28.3±0.6)E3	104.9±9.6	(80.2±1.6)E3	172.6±9.8	248.9±3.1	575.5±16.8 1
	ϵ abs = 10 -6

TABLE IV :

 IV Dense QP benchmarks (average runtime per time-step (IK) and total simulation runtimes).

	Task name	PROXQP	QUADPROG	OSQP	QPOASES	SCS	QPSWIFT	MOSEK
	UR3 IK	13.2±0.1µs	17.0±3.3µs	15.8±0.2µs	21.7±0.3µs	✗ 1	52.2±1.2µs	323.1±74µs 2
	UR5 IK	11.9±0.1µs	16.8±0.2µs	16.5±0.2µs	21.4±0.2µs	✗ 1	53.8±0.8µs	310.8±5.3µs 2
	DUAL ARMS IK	17.2±0.1µs	23.3±0.2µs	40.4±0.6µs	330.4±5.8µs	81.5±0.6µs	152.1±1.1µs	554.4±25.1µs 2
	KINOVAGEN2 IK	15.8±0.1µs	18.9±0.2µs	17.0±0.2µs	31.4±0.4µs	46.5±1.3µs	53.7±0.4µs	375.4±14.9µs 2
	SIGMABAN IK	14.1±0.2µs	25.2±0.4µs	45.2±0.9µs	523.6±4.8µs	68.6±0.5µs	224.9±2.0µs	452.5±8.8µs 2
	STRETCH IK	26.9±0.3µs	36.7±0.6µs	53.1±0.9µs	212.6±0.8µs	455.3±23.8µs 2	152.8±1.5µs	✗ 3
	CHAIN80 SQP	15.6±0.3ms	355.8±0.9ms	456.5±2.6ms	182±2.9ms	837±16.0ms	2193.9±22.3ms	1554.6±19.5ms 2
	CHAIN80W SQP	265.7±2.9ms	610.1±5.7ms	2141.9±22.4ms 467.4±3.3ms 1	✗ 4	✗ 4	3444.2±30.3ms 2
	HUMANOID MPC ϵABS = 10 -2	1.6±0.01ms	4.5±1.8ms	1.8±0.01ms	18.0±5.3ms	433.0±21.7ms 5	✗ 3	70.7±0.7ms
	HUMANOID MPC ϵABS = 10 -4	2.6±0.02ms	4.5±1.8ms	2.7±0.03ms	18.0±5.3ms	80.2±2.9ms 5	✗ 3	68.3±0.3ms
	HUMANOID MPC ϵABS = 10 -6	4.4±0.05ms	4.5±1.8ms	4.4±0.05ms	18.0±5.3ms	64.3±2.4ms 5	✗ 3	69.7±0.8ms

TABLE V :

 V Humanoid locomotion MPC problems with perturbations. Deviation from the steady state, the velocities of all masses, and the control action are penalized via the objective function. It forms a sequence of 101 QPs with 240 variables with only box constraints over the input variables. The second test problem CHAIN80W considers the same settings by adding also state constraints into the optimization problem in order to ensure that the chain does not hit a vertical wall close to the steady state. It results thus in a sequence of 101

	Noise Level	PROXQP	QUADPROG	OSQP	QPOASES	SCS	QPSWIFT	MOSEK
	10.0	11.9±9.7%	1.0±0.2%	1.9±0.9%	1.0±0.2%	1.0±0.2%	0±0.0%	1.0±0.2%
	5.0	58.38±36.4%	1.1±0.3%	2.1±1.1%	1.1±0.3%	1.1±0.4%	0±0.0%	1.1±0.3%
	1.0	100±0.0%	1.4±0.8%	3.5±2.4%	1.4±0.8%	1.5±1.1%	0±0.0%	1.4±0.9%
	0.5	100±0.0%	1.8±1.2%	5.5±3.8%	1.9±1.5%	2.1±1.6%	0±0.0%	1.8±1.2%
	0.1	100±0.0%	3.3±2.6%	51.6±36.7%	4.3±3.8%	4.9±4.3%	0±0.0%	3.3±2.6%
	0.05	100±0.0%	3.5±3.2%	97.6±13.5%	5.0±6.9%	7.7±6.5%	0±0.0%	3.5±3.2
	0.01	100±0.0%	4.4±4.4%	100±0.0%	7.7±9.8%	60.2±37.8%	0±0.0%	4.4±4.5%
	10 -3	100±0.0%	5.0±5.2%	100±0.0%	11.4±12.5%	100±0.0%	0±0.0%	5.0±5.2%
	10 -4	100±0.0%	5.0±5.2%	100±0.0%	15.5±16.8%	100±0.0%	0±0.0%	5.0±5.2%
	10 -5	100±0.0%	5.0±5.2%	100±0.0%	83.0±36.5%	99.1±8.9%	0±0.0%	5.1±5.3%
	10 -7	100±0.0%	5.0±5.2%	100±0.0%	100±0.0%	97±14.8%	0±0.0%	44.8±34.2%
	10 -9	100±0.0%	5.0±5.2%	100±0.0%	100±0.0%	100±0.0%	0±0.0%	100±0.0%
	0.0	100±0.0%	100±0.0%	100±0.0%	100±0.0%	100±0.0%	0±0.0%	100±0.0%
	chain is fixed on a wall while the three velocity components
	of the other end are used as control input with fixed lower and
	upper bounds. The prediction horizon of 16 seconds is divided
	into 80 control intervals. The model equations are derived
	from the linearisation of the non-linear ODE model (with 57
	states) at the steady state.

  , PROXQP solves the entire trajectory for noise amplitudes lower than 1.0 m/s 2 , while the second best solver OSQP reaches this level of performance for a noise amplitude of 10 -2 m/s 2 . All other solvers that are based on interior-point or active-set methods exhibit a less robust behavior on this task.

Real-robot balancing. We test PROXQP on real hardware with closed-loop model predictive control experiments on the

TABLE VI :

 VI Shifted geometric means (SGM) and failure rates (FR) for solving sparse QPs. The lower the better.

			PROXQP	QUADPROG	OSQP	QPOASES	SCS	QPSWIFT	MOSEK
		SGM		1		18.8	2.6	36.7	9.0	261.2	3279.8
		FR		0.0%		0.02%	0.0%	0.0%	0.02%	3.31%	25.38%
		1.0
	Ratio of problems solved	0.2 0.4 0.6 0.8				PROXQP OSQP SCS qpswift MOSEK quadprog qpOASES
		10 0 0.0	10 1	10 2 Performance ratio	10 3	10 4
	Fig.

Inria -Département d'Informatique de l' École normale supérieure, PSL Research University.

École des Ponts, Marne-la-Vallée, France.

ACKNOWLEDGMENTS

The authors would like to thank Viviane Ledoux for participating in the real-robot experiments. This work was supported in part by the French government under the management of Agence Nationale de la Recherche through the INEXACT project (ANR-22-CE33-0007-01) and as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). This work was also supported by the European Union through the AGIMUS project (GA no.101070165). Views and opinions expressed are those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.


Mélanie Munch
email: melanie.munch@gmail.com

Patrice Buche
email: patrice.buche@inrae.fr

Cristina Manfredotti
email: cristina.manfredotti@agroparistech.fr

Liliana Ibanescu
email: liliana.ibanescu@agroparistech.fr

An ontology-based pipeline to support the design of technical itineraries: application to composite food packaging

Keywords: Ontology, Knowledge Base, Probabilistic model, Causality, Reverse engineering

Due to the complexity of their data collection and their multi-scale approach, technical itineraries are complex and challenging to model. In this short abstract, we present a complete pipeline that aims to answer two kinds of expert competency questions: some by directly querying the data, and the others by learning a model able to reason with this complexity. By answering those, POND allows us to discover new applications and reverse engineer possible modifications.

POND: PO² ONtology Discovery

Composite food packaging is a challenging technical itinerary to design, due to the multiple parameters to be processed and optimized. Because of its potential variability, strong methods are needed to handle both its uncertainty and the need for structuration for its parameters.

In this short abstract, we present the POND (Process and Observation ONtology Discovery) pipeline, focused on the integration and the reasoning on knowledge from different sources. Based on the PO 2 Ontology [START_REF] Ibanescu | Po2-a process and observation ontology in food science. application to dairy gels[END_REF] (dedicated to represent Transformation Processes), it allows the representation and the probabilistic reasoning on complex technical itineraries.

POND is composed of three different steps 1:

• Step 1. Knowledge Collection. Knowledge comes initially from experimental data and experts interviews. Most of it is integrated to the PO² knowledge graph, which allows to semantize and structure complex and heterogeneous data from diverse sources. This structuration is exported as an RDF Graph that can then be queried through the SPARQL language. The expert has then the possibility of expressing Competency Questions (CQ) [START_REF] Uschold | Ontologies: principles, methods and applications[END_REF], which are natural-language questions that outline the scope of knowledge represented by an ontology. In our case, CQs can either be answered using a SPARQL query, or through a probabilistic reasoning with a specific model whose learning is detailed in the following steps. • Step 2. Mapping between PO² and PRM. Some CQs can only be answered through a probabilistic reasoning. This requires the introduction of the Probabilistic Relational Model (PRM) [START_REF] Koller | Probabilistic Graphical Models: Principles and Techniques -Adaptive Computation and Machine Learning[END_REF], which easily allows knowledge integration while robustly modeling uncertainty. The key idea is to exploit the expert knowledge defined in Step 1 to narrow the search space and guide the learning towards a model close to the reality in order to allow probabilistic and causal reasoning. To achieve the integration of this expert knowledge, two sub-steps have been defined: the Attributes Mapping where the expert defines the ontology's attributes relevant for the problem; and the Precedence Constraints Definition where the expert gives potential causal relations that will influence the learning [START_REF] Munch | Towards interactive causal relation discovery driven by an ontology[END_REF]. All of the information is then encompassed in a Relation Schema (RS), which is used to guide the PRM learning [START_REF]Introduction to Statistical Relational Learning (Adaptive Computation and Machine Learning)[END_REF]. • Step 3. PRM Exploitation for Reasoning. Once the PRM learned, it then becomes possible to analyse the learned model in order to answer the expert's CQ.

Application

Composite food packaging's goal is to provide new formulations of packages through the mix of bio-masses and polymeric matrices. In our application, a potential CQ is the possibility of reverse-engineer possible bio-masses that could be used to enhance the final package's features. In this case, the expert expresses constraints they would like to apply to the final products (such as better mechanical properties, for instance); from the PRM, it is then possible to deduce, through inference, scenarios (i.e. a set of conditions on the inputs) that would enhance the odds of satisfying the expert's constraints. These scenarios usually require a particular composition of the biomass, which can be used to query through different databases in order to discover and recommend new material that haven't been tested in the experiment yet.

Conclusion and Perspectives

In conclusion, we have in this short abstract presented a complete ontology-based pipeline that can integrate expert's knowledge in order to answer different kind of competency questions through two different means: 1) by SPARQL querying and 2) by learning a probabilistic model. This last solution offers a brand new look on the data, and even suggests other designs for the studied technical itinerary. On future works, it could be interesting to see how these discoveries can help formulate explanatory rules for the studied domain.

Figure 1 :

 1 Figure 1: Presentation of the POND Pipeline, which combines ontologies for data representation and structuration, and Probabilistic Relational Models for reasoning under uncertainty.