Nabiil Auckburally

Souheila El Bir

David Evrard
email: david.evrard@univ-tlse3.fr

Brigitte Dustou

Laure Latapie

Laurent Boumati

Christine Roques

Pierre Gros

Elucidation of the stoichiometry of Fe(III) and Fe(II) complexes formed with 2-(5-bromo-2-pyridylazo)-5-diethylaminophenol (5-Br-PADAP). Application to Fe(III) and Fe(II) Electrochemical Determination at Trace Level in Natural Water

Keywords: Complex stoichiometry, 5-Br-PADAP complexation, iron speciation, electrochemical sensor, natural water samples

This work provides some new information on the complexes formed between Fe(III) and Fe(II) with 2-(5-bromo-2-pyridylazo)-5-diethylaminophenol (5-Br-PADAP). The electrochemical behavior of both complexes was studied and their stoichiometry, which was uncertain till now, was definitely clarified as a 1:1 metal:ligand ratio on the basis of both electrochemical and spectrophotochemical data. The stability constants  for both complexes were also determined and found to be in the range 10 5 and 10 4 for Fe(III)/ and Fe(II)/5-Br-PADAP complex, respectively. Using 5-Br-PADAP as a complexing agent, the determination of Fe at trace level in water was achieved by adsorptive differential pulse voltammetry in the range 1 -10 µM. Under optimized conditions, a sensitivity of 0.0376 and 0.0063 µA µM -1 was found for Fe(III) and Fe(II), respectively. Tests on natural water from Canal du Midi were performed, and allowed values of 2.3 µM for Fe(III) and 1.7 µM for Fe(II) to be obtained. The corresponding value for total Fe was slightly lower than that found with ICP-OES measurements (ca. 4.0 vs. 4.5 µM, respectively), suggesting that 5-Br-PADAP is not a complexing agent strong enough to catch Fe cations in the complexes formed with organic matter nor to complex quantitatively free Fe(III) or Fe(II) species.


domain: phys spi

text:
Jonas Hamperl

Jean-Baptiste Dherbecourt
email: jean-baptiste.dherbecourt@onera.fr

Julien Totems

Martin Kjell

Carlota Mølster

Cherrie Canalias

Andrius Lee

Valdas Zukauskas

Rosa Pasiskevicius

Jean-Michel Santagata

Melkonian

Kjell Martin Mølster

Carlota Canalias

Cherrie Lee

Andrius Zukauskas

Valdas Pasiskevicius

Rosa Santagata

Jean-Michel Melkonian

Antoine Godard

Patrick Chazette

Cyrille Flamant

Myriam Raybaut

Preliminary range

Preliminary range-resolved detection of stable water isotopologues by differential absorption lidar using a 2 µm parametric source 1. Introduction

Differential absorption lidar (DIAL) is a remote gas sensing technique based on the use of tunable narrow-linewidth sources emitting radiation in and out of tune with the absorption line of the species of interest. To perform sensitive, range-resolved gas measurements, these sources need to deliver high peak power in the nanosecond regime [START_REF] Orr | Infrared LIDAR Applications in Atmospheric Monitoring[END_REF]. A potential application for such a device is the measurement of stable water isotopologues in the lower atmosphere, which would lead to improved understanding of processes governing the global hydrological cycle, with the prospects of improving the accuracy in the predictions made by atmospheric general circulation models [START_REF] Steen-Larsen | Evaluating the skills of isotope-enabled general circulation models against in situ atmospheric water vapor isotope observations[END_REF]. So far, despite the latest advances on multispecies DIALs [START_REF] Wagner | Multi-frequency differential absorption LIDAR system for remote sensing of CO2 and H2O near 1.6 μm[END_REF] and a first, recently reported demonstration of the possibility to probe isotopic water vapor lines [START_REF] Yu | Photon-counting distributed free-space spectroscopy[END_REF], no DIAL instrument has been able to quantify the relative abundance in HD 16 O (D), which requires collocated measurements of both isotopologues H2 16 O and HD 16 O. This can however be realized around 1.98 µm, as we theoretically demonstrated in [START_REF] Hamperl | Differential absorption lidar for water vapor isotopologues in the 1.98 µm spectral region: sensitivity analysis with respect to regional atmospheric variability[END_REF]. Based on these considerations, we realized a parametric source operating at 1.98 µm and performed a first demonstration of range-resolved probing of the water isotopologues H2 16 O and HD 16 O.

Parametric source and DIAL setup

The laboratory parametric source architecture is pictured in Fig. 1a. It is based on a nested-cavity optical parametric oscillator (NesCOPO) [START_REF] Godard | Nested Cavity Optical Parametric Oscillators -A Tunable Frequency Synthesizer for Gas Sensing[END_REF] combined with optical parametric amplification (OPA) stages, pumped by a commercial Nd:YAG laser (Innolas), delivering 12 ns-pulses at a repetition rate of 150 Hz. The NesCOPO allows the emission of a single-frequency radiation, tunable in a 1 nm-window around 1982.5 nm suitable for addressing H2 16 O and HD 16 O. The NesCOPO emitted wavelength is pulse-to-pulse modulated and recorded with a High Finesse WS6 IR wavemeter (Fig. 2). For each targeted species (H2 16 O and HD 16 O), such pulse-to-pulse wavelength modulation is necessary to emit the two wavelengths (namely ON and OFF) in and out of coincidence with the absorption line. This is performed by adjusting the NesCOPO cavity length using piezoelectric transducers. The OPA stage is based on three high-aperture (5x5 mm²) PPKTP crystals [START_REF] Hamperl | High Energy Parametric Laser Source and Frequency-Comb-Based Wavelength Reference for CO2 and Water Vapor DIAL in the 2 µm Region: Design and Pre-Development Experimentations[END_REF]. With a pump energy of 70 mJ, and idler rejection between the last two crystals, it was possible to obtain output signal energies in the range of 6-9 mJ. A transverse profile of the output beam after passing through a beam expander is shown in Fig. 1b. It is Gaussian with a measured M 2 <1.5 and a divergence full-angle of 300 µrad. The lidar receiver consists of a 10-inch, custom-made Newtonian telescope which focuses the backscattered radiation onto a 900 µm-multimode fiber, whose output is imaged onto a commercial InGaAs amplified photodiode (300 µm diameter). The telescope full-field-of view angle is 1.2 mrad. This DIAL experimental demonstration on water isotopologues is a first important step towards range-resolved D quantification in the vertical direction. For this purpose, future work will focus on the characterization of the DIAL errors and biases as well as on increasing the power performance of the transceiver, as outlined in [START_REF] Hamperl | High Energy Parametric Laser Source and Frequency-Comb-Based Wavelength Reference for CO2 and Water Vapor DIAL in the 2 µm Region: Design and Pre-Development Experimentations[END_REF], for higher sensitivity.

Fig. 1 .Fig. 2 .

 12 Fig. 1. (a) Experimental bench of the parametric laser source. (b) Transverse beam profile at the laser bench output

3 .Fig. 3 .

 33 Fig. 3. Raw (no energy normalization) averaged DIAL signals over 30 min and retrieved VMR for H2 16 O (a) and HD 16 O (b).

© 2 0 2 2 T h e Au t h o r ( s )

Aknowledgments

This research has received funding from the EU H2020 research and innovation program under grant agreement Nr. 821868, and from the French Agence Nationale pour la Recherche (ANR) under grant Nr. ANR-16-CE01-0009.


Oumayma Bounou

Jean Ponce

Justin Carpentier

Leveraging Proximal Optimization for Differentiating Optimal Control Solvers

Over the past few years, differentiable optimization has gained interest within machine learning, control, and robotics communities. It consists in computing the derivatives of the solutions of a given optimization problem which can then be used in learning algorithms. Until now, dedicated approaches have been proposed to compute the derivatives of various optimization problems (LPs, QPs, SOCPs, etc.). However, these approaches assume in general well-conditioned problems, limiting de facto their application to general optimal control problems (OCPs) widely used in robotics. In this work, we focus on the differentiation of solutions to such problems. We notably introduce a differentiable proximal formulation of equality-constrained LQR problems that accurately solves rank-deficient problems. This allows us to compute accurate gradients even in the case of problems that do not satisfy the standard linear independence constraint qualification (LICQ). Because any optimal control problem can be cast as an equalityconstrained LQR problem in the vicinity of the optimal solution, we show that our robust LQR derivative computation can be exploited to obtain the derivatives of general optimal control problems. We demonstrate the effectiveness of our approach in dynamics learning and parameter identification experiments in both linear and nonlinear optimal control problems.

I. INTRODUCTION

Trajectory optimization is a key part of programming robot motions. The task is encoded as a cost term, and physical constraints are encoded as path constraints. Over the past decades, differential dynamic programming (DDP) [START_REF] Mayne | Differential Dynamic Programming-A Unified Approach to the Optimization of Dynamic Systems[END_REF] and the iterative linear quadratic regulator (iLQR) [START_REF] Li | Iterative linear quadratic regulator design for nonlinear biological movement systems[END_REF] have become widespread and tractable approaches for solving complex robotic problems, ranging from UAVs navigation [START_REF] Marti-Saumell | Squash-box feasibility driven differential dynamic programming[END_REF] to legged locomotion [START_REF] Budhiraja | Differential dynamic programming for multi-phase rigid contact dynamics[END_REF]. While DDP was originally designed for unconstrained problems, variants have been proposed to account for diverse types of path constraints, ranging from simple bounds on the control inputs [START_REF] Tassa | Control-limited differential dynamic programming[END_REF] to equality [START_REF] Giftthaler | A projection approach to equality constrained iterative linear quadratic optimal control[END_REF], [START_REF] Kazdadi | Equality constrained differential dynamic programming[END_REF] and inequality constraints [START_REF] Lantoine | A hybrid differential dynamic programming algorithm for constrained optimal control problems. part 1: Theory[END_REF], [START_REF] Howell | Altro: A fast solver for constrained trajectory optimization[END_REF], [START_REF] Jallet | Constrained Differential Dynamic Programming: A primal-dual augmented Lagrangian approach[END_REF].

With recent progress in the development of advanced computational frameworks in machine learning (e.g., PyTorch [START_REF] Paszke | Pytorch: An imperative style, highperformance deep learning library[END_REF], TensorFlow [START_REF] Abadi | TensorFlow: Largescale machine learning on heterogeneous systems[END_REF], JAX [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF]), differentiable optimization has emerged as a generic approach for computing the derivatives of computational layers involved in mathematical programming problems. It can be used, for instance, to find the optimal design parameters of a robotic mechanism (aka codesign) given a target task [START_REF] Fadini | Computational design of energy-efficient legged robots: Optimizing for size and actuators[END_REF] or to estimate the sensitivity of an optimal solution to the parameters of the problem [START_REF] Oshin | Parameterized Differential Dynamic Programming[END_REF]. Differentiable optimization generally relies on differentiating the optimality conditions (e.g., the Karush-Kuhn-Tucker 1 Inria 2 Inria, Département d'Informatique de l'École normale supérieure, (ENS-PSL, CNRS, Inria) 3 Courant Institute and Center for Data Science, New York University conditions) associated with a given problem. In the field of convex optimization, effective approaches are now available for computing derivatives of solutions to standard quadratic programming (QP) instances [START_REF] Amos | Optnet: Differentiable optimization as a layer in neural networks[END_REF], LQR problems [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], and more general convex programming problems [START_REF] Agrawal | Differentiable convex optimization layers[END_REF], which are very common in control and robotics. However, when deriving the KKT conditions, these approaches often assume that they are well-conditioned at the optimum (e.g., linearly independent constraints qualification, low condition number of the KKT matrix), which cannot be ensured for a large majority of control problems (e.g., inverse or forward dynamics of legged robots making redundant contacts with their environment or singular configurations of robotic arms). This is even more true in machine learning applications, where optimization stages are used as differentiable layers, and the intermediate problems that are being solved during the learning process are not guaranteed to satisfy the required optimality qualifications since they are never enforced during training.

We propose to overcome these limitations by leveraging the proximal method of multipliers [START_REF] Rockafellar | Augmented Lagrangians and applications of the proximal point algorithm in convex programming[END_REF] to evaluate the derivatives of solutions to optimal control problems robustly. In particular, we develop a generic solution to the computation of these derivatives in the case of the equality-constrained linear quadratic regulator (LQR), which is, as highlighted in [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], the core block required to compute the derivatives of solutions to general nonlinear optimal control problems. Following the experimental validation procedure of [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], we demonstrate that unlike existing solutions, our proximal approach properly handles ill-conditioned system identification problems. We also show that our framework can be used on top of an external nonlinear optimal control solver, enabling our approach to operate on nonlinear problems. This paper is organized as follows. We first introduce our proximal formulation of equality-constrained LQR problems (Sec. II). We then present our proximal method for differentiating through these problems (Sec. III), and the extension of our approach to the nonlinear case (Sec. IV). In section V, we evaluate and benchmark the proposed solution against alternative solutions of the state of the art on various system identification problems of the literature.

II. PROBLEM STATEMENT AND PROXIMAL RESOLUTION OF EQUALITY-CONSTRAINED LQR

A. Problem statement

We address the equality-constrained LQR problem of finding the optimal sequence of states and controls that solve

min X,U T -1 t=0 l t (x t , u t ) + l T (x T ), s.t    x t+1 = f t (x t , u t ) x 0 = x * 0 c t (x t , u t ) = 0 , (1)
where

X = x T 0 ... x T T T , U = u T 0 ... u T T -1
T , l t is a quadratic cost function, f t is a linear dynamics function and c t are path constraints. We define the running cost:

l t (x t , u t ) = 1 2 (x T t Q t x t + u T t R t u t ) + q T t x t + w T t u t , (2)
the terminal cost:

l T (x T ) = 1 2 (x T T Q f x T ) + q T T x T , (3)
the dynamics:

f t (x t , u t ) = A t x t + B t u t + d t , (4)
and the path constraints

C t x t + D t u t -e t = 0 and x 0 = x * 0 . (5)
In the case of LQR problems, several approaches [START_REF] Boyd | Lqr[END_REF] can be used to solve problem [START_REF] Mayne | Differential Dynamic Programming-A Unified Approach to the Optimization of Dynamic Systems[END_REF]. We follow here the Bellman principle of optimality [START_REF] Liberzon | Calculus of variations and optimal control theory: a concise introduction[END_REF]. Starting from t = T , we introduce the value function V t computed recursively backward in time as:

V t (x t ) = min ut,xt+1 l t (x t , u t ) + V t+1 (x t+1 ) (6)

B. Multiple-shooting

Following [START_REF] Jallet | Implicit Differential Dynamic Programming[END_REF], we also relax the dynamics constraints under Eq. ( 6) by using multiple shooting [START_REF]A multiple shooting algorithm for direct solution of optimal control problems[END_REF], which consists in adding auxiliary state variables in order to stabilize the optimization procedure, especially for numerically or intrinsically unstable systems. We introduce an auxiliary variable y t , such that Eq. ( 6) maps to:

V t (x t ) = min ut,yt l t (x t , u t ) + V t+1 (y t ) (7)
s.t. y t = f t (x t , u t ) and c t (x t , u t ) = 0.

The new problem converges to the solution of problem [START_REF] Giftthaler | A projection approach to equality constrained iterative linear quadratic optimal control[END_REF] as demonstrated by Schmidt et al. in [START_REF] Schmidt | Convergence rates of inexact proximal-gradient methods for convex optimization[END_REF]. In the sequel, we use the usual shorthands V x and V xx to denote respectively the first and second derivatives of V t , and V ′ to denote V t+1 .

At time t, we solve problem [START_REF] Kazdadi | Equality constrained differential dynamic programming[END_REF] by solving the equivalent min max problem:

V t (x t ) = min yt,ut max λt,νt L t (x t , u t , y t , λ t , ν t ), (8)
where L t is the Lagrangian of the problem:

L t (x t , u t , y t , λ t , ν t ) = l t (x t , u t ) + V ′ (y t ) + λ T t (y t -f t (x t , u t )) (9)
+ ν T t c t (x t , u t ). Here, y t and u t are the primal variables and λ t and ν t are the dual ones. In the LQR setting, V T (x T ) = 1 2 x T T Q f x T + q T T x T , hence V t is also quadratic for each t [START_REF] Liberzon | Calculus of variations and optimal control theory: a concise introduction[END_REF]. The LQR problem can then be solved using the dynamic programming principle by solving at each time t the linear system obtained by deriving the KKT conditions at optimality:

    R t 0 B T t D T t 0 V ′ xx -I 0 B t -I 0 0 D t 0 0 0         u t y t λ t ν t     = -     w t V ′ x A t x t + d t C t x t -e t     , ( 10
)
where the leftmost matrix is the so-called KKT matrix.

Because of the multiple shooting approach, this system has an additional variable y t . However, the KKT matrix has a blocksparse structure, which can be efficiently exploited when solving the system. It should be noted that such a matrix can have very poor conditioning, especially when V ′ xx or R t is not positive definite or when D t is rank-deficient. As mentioned in the introduction, such a scenario is likely to happen in practice, especially in a machine learning context where the cost matrices R t and Q t and the path constraints matrix D t are learnable parameters. To cope with these limitations, we propose solving a proximal reformulation of this problem, which we detail in the following.

C. Proximal reformulation

As shown by Rockafellar in [START_REF] Rockafellar | Augmented Lagrangians and applications of the proximal point algorithm in convex programming[END_REF], solving the min-max problem ( 8) is equivalent to solving the equivalent problem

min yt,ut max λt,νt L t,ρµ (x t , u t , y t , λ t , ν t ), (11)
where ρ and µ are positive constants,

L t,ρµ (x t , u t , y t , λ t ) = l t (x t , u t ) + V t+1 (y t ) + λ T t (y t -f (x t , u t )) + ν T t c t (x t , u t ) + ρ 2 ∥u t -u - t ∥ 2 2 + ρ 2 ∥y t -y - t ∥ 2 2 - µ 2 ∥λ t -λ - t ∥ 2 2 - µ 2 ∥ν t -ν - t ∥ 2 2 , ( 12
)
and the four L 2 regularization terms vanish over the iterations of a backward/forward procedure for solving problem [START_REF] Paszke | Pytorch: An imperative style, highperformance deep learning library[END_REF] for all t. The so-called proximal regularization [START_REF] Parikh | Proximal algorithms[END_REF] of the primal variables u t , y t and dual variables λ t , ν t is used for solving problems that are not both primally and dually strictly convex. It is central to our approach and will be notably used to solve constrained LQR problems and estimate the sensitivity of their solutions to problem parameters robustly.

Let us now present the backward/forward procedure associated with our approach. Backward pass. For all t backward in time, the new KKT conditions at time t are given by:

    R t + ρI 0 B T t D T t 0 V ′ xx + ρI -I 0 B t -I -µI 0 D t 0 0 -µI         du t dy t dλ t dν t     = -     R t u - t + B T t λ - t + D T t ν - t + w t V ′ xx y - t + V ′ x -λ - t A t x t + B t u - t + d t -y - t C t x t + D t u - t -e t     , (13)
where dv = v -v -, for v in {u t , y t , λ t , ν t }. When solving Eq. ( 13) backward in time, x t is unknown, but du t , dy t , dλ t and dν t can be expressed as affine functions of

x t :        du t = Γ t x t + γ t , dy t = M t x t + m t , dλ t = Ξ t x t + ξ t , dν t = Ω t x t + ω t , , (14)
and the coefficients Γ t , γ t , M t , m t , Ξ t and ξ t can be obtained by solving the linear system:

K t     Γ t γ t M t m t Ξ t ξ t Ω t ω t     = -     0 R t u - t + B T t λ - t + D T t ν - t + w t 0 V ′ xx y - t + V ′ x -λ - t A t Bu - t + d t -y - t -C t D t u - t -e t     , (15)
where K t is the KKT matrix from Eq. ( 13). In the backward pass of the algorithm, we compute the coefficients from Eq. ( 15) and update the value function's first and second derivatives V x and V xx for each t (which is straightforward since V t is quadratic).

Forward pass. We update u t , y t , λ t and ν t for all t, forward in time.

       u + t = u - t + du t , y + t = y - t + dy t , λ + t = λ - t + dλ t , ν + t = ν - t + dν t . . ( 16
)
We thus solve the LQR problem by iteratively repeating this backward/forward procedure until convergence. Our stopping criterion uses the infinity norm of the gradient of L ρ,µ

∇ y,u L ρµ (y, u) ∇ λ,ν L ρµ (λ, ν) ∞ ≤ ϵ, (17)
where L ρµ stacks L t,ρµ for all t and ϵ is a fixed tolerance parameter.

III. PROXIMAL DIFFERENTIATION OF EQUALITY-CONSTRAINED LQR Our objective is to compute derivatives of LQR outputs (i.e., the trajectories x and u) with respect to the problem parameters (i.e., the dynamics and cost matrices A t , B t , Q t , R t , the path constraints C t , D t , d t , e t , and the initial condition x * 0 ) to be able to plug the LQR solver as a differentiable layer in a trainable model. Following [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], we first reformulate [START_REF] Giftthaler | A projection approach to equality constrained iterative linear quadratic optimal control[END_REF] as the minimization of an equivalent QP problem of the form:

min X 1 2 X T HX + q T X, s.t ÂX = b, (18)
with the following augmented matrices:

H = diag(Q t , ..., R t , ...), q = [q T 0 ... q T T w T 0 ... w T T -1 ] T , Â =               -I 0 0 ... 0 0 0 ... 0 A 0 -I 0 ... 0 B 0 0 ... 0 0 A 1 -I ... 0 0 B 1 ... 0
0 D T -1 ... 0 0 0 ... ... 0 C T 0 ... 0 0               , X = [x T 0 ... x T T u T 0 ... u T T -1 ] T , and b = -[x * T 0 d T 1 ... d T T e T 0 ... e T T ] T .
The KKT conditions associated with problem 18 at the optimum correspond to:

K X * Λ * = -q b , where K = H ÂT Â 0 (19)
and Λ * is a vector of multipliers associated with the constraints ÂX = b and containing the stacked multipliers λ t and ν t from Eq. ( 9). for all t. We drop the " * " superscript for clarity in the following.

In machine learning problems, we are looking for derivatives of the form ∂r/∂p, where r is a scalar function and p is a parameter of the problem . Here, p lies in P t = {A t , B t , Q t , q t , w t , C t , D t , x * 0 , d t , e t }. To obtain ∂r/∂p directly, we use the same trick as in [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF] and [START_REF] Le Lidec | Differentiable simulation for physical system identification[END_REF], and obtain expressions for the partial derivatives

∂r/∂At = λt+1G T x,t + xtG T λ,t+1 , ∂r/∂Bt = λt+1G T u,t + G λ,t+1 u T t , ∂r/∂Qt = xtG T x,t , ∂r/∂Rt = Gu,tu T t , ∂r/∂Ct = -xtG T ν,t , ∂r/∂et = Gν,t, ∂r/∂Dt = -utG T ν,t , ∂r/∂dt = G λ,t , ∂r/∂qt = Gx,t, ∂r/∂wt = Gu,t, (20)
as functions of a vector G that verifies

KG = Z, ( 21
)
where Z is the vector obtained by stacking all the vectors ∂r/∂v t with v t in {x t , u t , λ t , ν t }, for all t. The quantities in Z are assumed to be known. In fact, they correspond to the derivatives of r with respect to the outputs of the LQR, i.e., the derivatives of the computational block following the LQR block in the considered computational graph, and can be obtained with back-propagation. We can also write G by stacking the vectors G v,t for all t, with v t in {x t , u t , λ t , ν t } where each vector G v,t is the same size as v t . In [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], Amos et al. notice that G verifies the same equation that the vector [X T Λ T ] T verifies in Eq. [START_REF] Rockafellar | Augmented Lagrangians and applications of the proximal point algorithm in convex programming[END_REF]. Thus, G is the solution to a new LQR problem, and as such can be obtained efficiently without explicitly inverting the large KKT matrix K in Eq. [START_REF] Liberzon | Calculus of variations and optimal control theory: a concise introduction[END_REF]. Instead, we solve a similar LQR problem to the one to which x and u are solutions. The KKT matrices of this new problem, which we refer to as the LQR derivatives problem, are the same as the original ones. Only the values of the right-hand side term in the large QP formulation are modified, which correspond to the parameters q t , w t , d t , and e t in the LQR formulation. Their values for the LQR derivatives problem are:

q t = ∂r/∂x t w t = ∂r/∂u t d t = ∂r/∂λ t e t = ∂r/∂ν t . (22)
In [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], the authors assume the LICQ conditions are satisfied, whether it is to solve the LQR problem or to compute its derivatives, which implies that K t is always invertible. However, as discussed in the introduction, this limits the type of control problems the method can be applied to since it is not the case for many of them. Thus, unlike [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], we cope with these limitations by leveraging the proximal formulation introduced in section II-C.

We solve the LQR problem of which G x,t , G u,t , G λ,t and G ν,t are the solutions and multipliers by solving iteratively until convergence the following linear system for all t, backward in time starting from t = T :

K t     dG u,t dG y,t dG λ,t dG ν,t     = -     R t G - u,t + B T t G - λ,t + D T t G - ν,t + ∂r/∂u t V ′ xx G - y,t + V ′ x -G - λ,t A t G x,t + B t G - u,t + ∂r/∂λ t -G - y,t C t G x,t + D t G - u,t -∂r/∂ν t     , (23)
where

K t =    Rt + ρI 0 B T t D T t 0 V ′ xx + ρI -I 0 Bt -I -µI 0 Dt 0 0 -µI    (24)
and obtain new updates

    G u,t G y,t G λ,t G ν,t     =     G - u,t + dG u,t G - y,t + dG y,t G - λ,t + dG λ,t G - ν,t + dG ν,t     . (25)
This allows us to find the optimal variable G accurately, even in ill-conditioned problems (e.g., rank-deficient matrices K t ). This is extremely important in practice since the gradients of interest, ∂r/∂p where p is in P, are functions of the values of G (equation ( 20)). Thus, any inaccurate solution G would lead to back-propagating incorrect gradients ∂r/∂p, resulting in unstable training procedures.

IV. EXTENSION TO NONLINEAR CASES

When f , l, and c are nonlinear, problem (1) corresponds to an equality-constrained nonlinear optimal control problem. A standard method for solving it is the so-called iterative LQR (iLQR) algorithm [START_REF] Li | Iterative linear quadratic regulator design for nonlinear biological movement systems[END_REF], performing a cascade of LQR problems by linearizing the dynamics and making a quadratic approximation of the cost function around the nominal state and control trajectories. In the case of equality-constrained nonlinear optimal control problems, the approximations can be written as follows: iLQR-based approaches for solving unconstrained and constrained problems of this kind [START_REF] Tassa | Control-limited differential dynamic programming[END_REF], [START_REF] Howell | Altro: A fast solver for constrained trajectory optimization[END_REF], [START_REF] Jallet | Implicit Differential Dynamic Programming[END_REF] proceed iteratively until a convergence criterion is met. At the optimum, the KKT conditions of the nonlinear constrained optimal control problem (1) correspond to those of the last LQR problem solved in the iLQR. Thus, the derivatives of problem (1) in the general case can be obtained with the method described in III. The solution extends to OCPs containing inequality constraints, since active inequality constraints at the optimum play the role of pure equality constraints.

min δX,δU T -1 t=0 1 2 δx T t L t,
Finally, if we assume that the dynamics f , constraints c, and cost function l are parameterized by some parameter θ, then it follows that the set of first and second-order derivatives D l,f , and specifically the one at optimality D * l,f will also be parameterized by this parameter θ. In a learning framework, the derivative of interest, i.e., ∂r/∂θ, with r some scalar loss function over the parameters can thus be obtained by applying the chain rule, leading to:

∂r ∂θ = p∈D * l,f ∂r ∂p ∂p ∂θ , ( 28
)
where the derivatives ∂r/∂p can be obtained with the method described in paragraph III, and the derivatives ∂p/∂θ are obtained with backpropagation. Overall, our approach is quite general and we can use any off-the-shelf iLQR solver to get the set of derivatives of the OCP at optimality D * l,f to formulate the quadratic cost problem with linear equality constraints at optimality, which we then solve with our approach.

V. EXPERIMENTS

We first evaluate our method for solving LQR problems with terminal constraints. We then evaluate it on cost and parameter identification problems for both linear (LQR) and nonlinear systems (pendulum and cartpole). The code associated with this paper, written in Python, will be released as open-source.

A. LQR problems with terminal constraints

We solve an LQR problem with a terminal equality constraint (Sec. II-A) following the approach described in section II-C. We compare the solutions of both our solver and CVXPY [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF], [START_REF] Agrawal | A rewriting system for convex optimization problems[END_REF], which is based on OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]. We run both solvers on three sets of parameter sizes {n, d, T }, where Tab. I reports the average results of our experiments. In the first and second rows, our solver converges to a solution for almost every experiment, while CVXPY only converges to a solution for up to 32% of them on average. We explain this difference by the fact that for these problems, the KKT matrix involved in the last step of the dynamic programming problem has a condition number of 10 17 in the unregularized formulation, and of 10 7 in our regularized approach. In examples where both solvers converge, the proposed approach converges to a more accurate solution. This first set of results shows both increased robustness and better accuracy of the proximal solver against a more classic approach for solving LQR problems.

We also report the average results on well-conditioned cases (third row). In such cases, both solvers converge to equally good solutions, but CVXPY is much faster since our implementation is in an interpreted language (Python) that could easily be moved to a compiled one (e.g., C++) for much better efficiency.

B. System identification: LQR problems

We reproduce the system identification experimental setting of [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF]. Given optimal trajectories in states and controls of systems with linear dynamics and quadratic cost, the goal is to identify these dynamics and cost parameters. Formally, we solve the LQR problem of equations ( 2)-( 5) with q t = 0 and w t = 0: We observe M trajectories [x * ,i 0 , . . . x * ,i T ] and [u * ,i 0 , . . . , u * ,i T -1 ] (i in 1, . . . , M ) that are optimal solutions to problem [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF] with different initial conditions x * 0 (chosen at random). The states x t are vectors of size n (between 2 and 10), the controls u t are vectors of size d (between 3 and 10), and the control horizon is T (between 5 and 20). Our goal is to identify the dynamics matrices (A, B) and the cost matrices

min X,U 1 2 T -1 t=0 (x T t Qxt + u T t Rut) + 1 2 x T T Q f xT , s.t xt+1 = Axt + But for t in 0, . . . , N -1, x0 = x * 0 . (29)
(Q, Q f , R).
In other words, we want to estimate θ, where θ can be any of A, B, Q, Q f , R, or a combination of two or more of these matrices. We solve:

min θ M i=1 T t=0 ∥x * ,i t -x i t (θ)∥ 2 2 + ∥u * ,i t -u i t (θ)∥ 2 2 ( 30
)
where x i t and u i t are solutions of the LQR problem parameterized by θ with initial condition x * ,i 0 . Using the same notations as in [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], we define the optimal trajectory vector i, τ * ,i as

τ * ,i = x * ,i 0 T . . . x * ,i T T u * ,i 0 T . . . u * ,i T -1 T T . (31)
Here, τ * ,i is a vector of size p = (T +1)n+T d. The problem is now reduced to solving

min θ M p k=0 ∥r k (θ)∥ 2 2 , (32)
where r k (θ) = Γ * k -Γ k , τ * ,i k (respectively τ i k (θ)) is the k-th component of vector τ * ,i (respectively τ i ), and Γ * k (respectively Γ k ) contains stacked vectors τ * ,i (respectively τ i ). Problem [START_REF] Hinton | Rmsprop[END_REF] is a non-linear least-squares problem that can be solved using methods such as Gauss-Newton (GN) [START_REF] Floudas | Encyclopedia of Optimization[END_REF] or Levenberg-Marquardt (LM) [START_REF] Moré | Computing a trust region step[END_REF]. When a step in a solution of an optimal control problem boils down to least squares, stochastic gradient descent should not be used for this step (as expected from the optimization literature), which partly explains what happens in [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], where the identification experiments fail in half the trials when using gradient descent to optimize Eq. [START_REF] Floudas | Encyclopedia of Optimization[END_REF].

Table II shows results averaged on 100 experiments. An experiment is considered successful if the identification error ∥θ -θ * ∥ ∞ reaches the threshold 5.10 -6 in less than 50 iterations of the LM algorithm or 2000 epochs of RMSProp [START_REF] Hinton | Rmsprop[END_REF]. We see that the identification succeeds 85% of the time when using a least-squares method, while it only succeeds 20% of the time when using RMSProp, which demonstrates the ineffectiveness of stochastic gradient methods in this setting. Note that for the experiments using RMSProp, we have selected the hyperparameters (learning rate and batch size), which achieved the lowest errors. It should be noted, however, that implementing Levenberg-Marquardt-like methods requires computing the full Jacobian matrix of the residual function r (instead of just Jacobianvector products as in SGD methods), which is a function that scales linearly with T , n and d.

In the following experiments, the identification problem (32) was solved using the LM method for both our solver and the solver from [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF]. Figures 1 and2 show the system identification error and optimization loss as functions of the number of iterations using both the solver from [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF], which we refer to as diff-mpc, and our solver on 4 trials (100 experiments were run, each with a different set of parameters to identify, but we only show 4 randomly selected ones here). In this experiment, the parameters to identify are the dynamics and control matrices A and B. The control cost matrix R is set to the identity matrix, and the state cost matrix Q is set to 10 -4 I n . In this experiment, both solvers converge, but our regularized solver converges in half as many iterations on all the trials. Figures 3 and4 show the system identification error and optimization loss as functions of the number of iterations using both solvers again, on experiments where the parameter to identify is the matrix Q (the identity matrix). The approach of [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF] fails on all trials, demonstrating the importance of regularizing the LQR solver to avoid performing optimization steps with wrong gradients.

C. System identification: nonlinear dynamics

We demonstrate the effectiveness of our approach in identifying parameters of a nonlinear system from observed state and control trajectories. We consider a pendulum parameterized by its mass and its pole length and a cartpole parameterized by its pole mass, its cart mass, and its pole length. We seek to recover these parameters for both systems. The state of the pendulum at time t is defined by s t = [cos θ t , sin θ t , θt ], and the state of the cartpole by s t = [x t , ẋt , cos θ t , sin θ t , θt ] (translation over the x axis). We have conducted experiments with pendulums and cartpoles with 10 various sets of parameters and report the losses for 4 of them for each system (see table III). For the pendulum, the state cost matrix is set to Q = 10 -3 I n , for the cartpole, it is set to Q = 10 -1 I n . The control cost matrix is set to R = 10 -1 I d for both systems. A linear cost term is added to bring the pole upward and the cart centered at the origin. Here, the problem is also a nonlinear least-squares one, and it can be formulated as in Eq. ( 32), with θ = {m, l} for the pendulum, and θ = {m, m c , l} for the cartpole. We solve it using the LM method for both solvers and report the results in Figures 5 and6. We see that, unlike diff-mpc, our method converges to the ground truth parameters in all experiments. When both methods converge, ours converges in fewer iterations to a better solution.

VI. CONCLUSION

We have introduced a regularized differentiable equalityconstrained LQR solver with a generic formulation that handles path constraints. When used in learning frameworks, our solver is robust to ill-conditioned problems and shows better performance on system identification experiments compared to the unregularized approach introduced in [START_REF] Amos | Differentiable mpc for end-to-end planning and control[END_REF]. Since our formulation of the LQR problem is generic, it can also be applied to differentiate accurately through more general optimal control problems. Future work should focus on experiments on real robotic systems to demonstrate the effectiveness and robustness of this approach on real-world robotic tasks.

Fig. 1 :

 1 Fig. 1: Identification error on identifying the matrices A and B with Q = 10 -4 I n . Pairs of curves with the same colors are identification experiments on the same problem parameters solved using different solvers: diff-mpc in dashed lines and ours in solid lines. The same convention is used in all figures.

Fig. 2 :

 2 Fig. 2: Optimization loss on identifying the matrices A and B with Q = 10 -4 I n .

Fig. 3 :

 3 Fig. 3: Identification error. Identification of Q.

Fig. 4 :

 4 Fig. 4: Optimization loss. Identification of Q.

Fig. 7 :

 7 Fig. 7: Cartpole. Identification error. Identification of m, m c and l.

Fig. 8 :

 8 Fig. 8: Cartpole. Optimization loss. Identification of m, m c and l.

Fig. 5 :

 5 Fig. 5: Pendulum. Identification error. Identification of m and l.

Fig. 6 :

 6 Fig. 6: Pendulum. Optimization loss. Identification of m and l.

  s.t. x t+1 = f t (x t , u t ) and c t (x t , u t ) = 0, with the terminal condition V T (x T ) = l T (x T ). Here, Q t and R t are symmetric positive semi-definite matrices, A t , B t , C t and D t are general matrices, and e t and d t are vectors, for all t. x * 0 is the initial state. It is worth mentioning that C t and D t might be rank deficient.

  C t,x δx t + C t,u δu t + c 0 = 0, where δx t and δu t stand for the step directions and D l,f = {L t,xx , L t,uu , L t,xu , l t,x , l t,u , F t,x , F t,u , C t,x , C t,u } is the set of first and second-order derivatives of l, f and c respectively.

	xx δx t +	1 2	δu T t L t,uu δu t + x T t L t,xu δu t
	+ l T t,x δx t + l T t,u δu t			(26)
	s.t. F t,x δx t + F t,u δu t + f 0 = 0,	(27)
	s.t.

TABLE I :

 I Solvers performance. Comparison of LQR problems with terminal equality constraints. is the system dimension, d is the control dimension, and T is the time horizon. For each set of parameter sizes, we run 100 experiments with randomly generated time-invariant LQR problems. The dynamics matrix A is forced to have all singular values lower than 1, and B is a matrix with random coefficients sampled uniformly in [0, 1]. The cost matrices Q and R are set to respectively 10 -2 I n and 10 -1 I d . All experiments are run on a single CPU. A solver is considered successful when both primal and dual constraints are satisfied. Feasibility denotes the infinity norm on primal constraints (dynamics constraints ∥x * t+1 -Ax * t -Bu * t ∥ ∞ ). Distance to goal is the average infinite norm ∥x T -x * T ∥ ∞ , with x * T the target terminal constraint.

	Parameters	solver	success	feasibility distance to goal
	T = 20,	cvxpy	2	6.10 -6	6.10 -6
	n = 15, d = 3	ours	93	9.10 -9	9.10 -9
	T = 20,	cvxpy	32	2.10 -8	2.10 -8
	n = 10, d = 2	ours	98	4.10 -9	4.10 -9
	T = 20,	cvxpy	100	9.10 -12	9.10 -12
	n = 15, d = 5	ours	100	2.10 -14	3.10 -15

n

TABLE II :

 II Success rate of identification methods. We compare Levenberg-Marquardt and RMSProp optimization success rates (in %).

		T=5, n=3, d=3	T=20, n=3, d=3 T=10, n=5, d=2
	RMSProp	21	20	21
	LM	91	85	90