Neural Contraction Metrics for Robust Estimation and Control: A Convex Optimization Approach

This paper presents a new deep learning-based framework for robust nonlinear estimation and control using the concept of a Neural Contraction Metric (NCM). The NCM uses a deep long short-term memory recurrent neural network for a global approximation of an optimal contraction metric, the existence of which is a necessary and sufficient condition for exponential stability of nonlinear systems. The optimality stems from the fact that the contraction metrics sampled offline are the solutions of a convex optimization problem to minimize an upper bound of the steady-state Euclidean distance between perturbed and unperturbed system trajectories. We demonstrate how to exploit NCMs to design an online optimal estimator and controller for nonlinear systems with bounded disturbances utilizing their duality. The performance of our framework is illustrated through Lorenz oscillator state estimation and spacecraft optimal motion planning problems.


I. INTRODUCTION
P ROVABLY stable and optimal state estimation and control algorithms for a class of nonlinear dynamical systems with external disturbances are essential to develop autonomous robotic explorers operating remotely on land, in water, and in deep space. In these next generation missions, these robots are supposed to intelligently perform complex tasks with their limited computational resources, which are not necessarily powerful enough to run optimization algorithms in real-time.
Our main contribution is to introduce a Neural Contraction Metric (NCM), a global representation of optimal contraction metrics sampled offline by using a deep Long Short-Term Memory Recurrent Neural Network (LSTM-RNN) (see Fig. 1), and thereby propose a new framework for provably stable and optimal online estimation and control of nonlinear systems with bounded disturbances, which only requires one function evaluation at each time step. A deep LSTM-RNN [1], [2] is a recurrent neural network with an improved memory structure proposed to circumvent gradient vanishing [3] and is a universal approximator of continuous curves [4]. Contrary to previous works, the convex optimization-based sampling methodology in our framework allows us to obtain a large enough dataset of the optimal contraction metric without assuming any hypothesis function space. These sampled metrics, the existence of which is a necessary and sufficient condition for exponential convergence [5], can be approximated with  arbitrary accuracy due to the high representational power of the deep LSTM-RNN. We remark that this approach can be used with learned dynamics [6] as a nominal model is assumed to be given. Also, this is distinct from Lyapunov neural networks designed to estimate a largest safe region for deterministic systems [7], [8]: the NCM provides provably stable estimation and control policies, which have a duality in their differential dynamics and are optimal in terms of disturbance attenuation. The NCM construction is summarized as follows.
In the offline phase, we sample contraction metrics by solving an optimization problem with exponential stability constraints, the objective of which is to minimize an upper bound of the steady-state Euclidean distance between perturbed and unperturbed system trajectories. In this paper, we present a convex optimization problem equivalent to this problem, thereby exploiting the differential nature of contraction analysis that enables Linear Time-Varying (LTV) systems-type approaches to Lyapunov function construction. For the sake of practical use, the sampling methodology is reduced to a much simpler formulation than those of [9]- [11] derived for Itô stochastic nonlinear systems. These optimal contraction metrics are sampled using the computationally efficient numerical methods for convex programming [12]- [14] and then modeled by the deep LSTM-RNN as depicted in Fig. 1. In the online phase, contraction metrics at each time instant are computed by the NCM to obtain the optimal feedback estimation and control gain or a bounded error tube for robust motion planning [15], [16].
We illustrate how to design an optimal NCM-based estimator and controller for nonlinear systems with bounded disturbances, utilizing the estimation and control duality in differential dynamics analogous to the one of the Kalman filter and Linear Quadratic Regulator (LQR) in LTV systems. Their performance is demonstrated using Lorenz oscillator state estimation and spacecraft optimal motion planning problems.
Related Work: Contraction analysis, as well as Lyapunov theory, is one of the most powerful tools in analyzing the stability of nonlinear systems [5]. It studies the differential (virtual) dynamics for the sake of incremental stability by means of a contraction metric, the existence of which leads to a necessary and sufficient characterization of exponential stability of nonlinear systems. Finding an optimal contraction metric for general nonlinear systems is, however, almost as difficult as finding an optimal Lyapunov function.
Several numerical methods have been developed for finding contraction metrics in a given hypothesis function space. A natural application of this concept is to represent their candidate as a linear combination of some basis functions [17]- [20]. In [21], [22], a tractable framework to construct contraction metrics for dynamics with polynomial vector fields is proposed by relaxing the stability conditions to the sum of squares conditions. Although it is ideal to use a larger number of basis functions seeking for a more optimal solution, the downside of this approach is that the problem size grows exponentially with the number of variables and basis functions [23].
We could thus alternatively rely on numerical schemes for sampling data points of a Lyapunov function without assuming any hypothesis function space. This includes the state-dependent Riccati equation method [24], [25] and it is proposed in [9]- [11] that this framework can be improved to obtain an optimal contraction metric, which minimizes an upper bound of the steady-state mean squared tracking error for nonlinear stochastic systems. However, solving a nonlinear system of equations or an optimization problem at each time instant is not suitable for systems with limited online computational capacity. The NCM addresses this issue by approximating the sampled solutions by the LSTM-RNN.
II. PRELIMINARIES We use the notation x and A for the Euclidean and induced 2-norm, and A 0, A 0, A ≺ 0, and A 0 for positive definite, positive semi-definite, negative definite, and negative semi-definite matrices, respectively. Also, sym(A) = A + A T and I denotes the identity matrix. We first introduce some preliminaries that will be used to construct an NCM.

A. Contraction Analysis for Incremental Stability
Consider the following perturbed nonlinear system: where t ∈ R ≥0 , x : R ≥0 → R n , f : R n × R ≥0 → R n , and d : Theorem 1: Let x 1 (t) and x 2 (t) be the solution of (1) with d(t) = 0 and d(t) = 0, respectively. Suppose there exist Then the smallest path integral is exponentially bounded, thereby yielding a bounded tube of states: Rewriting this inequality using the relations √ ωR(t) ≥ x 2 x 1 δ x and 1 ≤ ω/ω ≤ ω/ω due to 0 < ω ≤ ω completes the proof. Input-to-state stability and finite-gain L p stability follow from (4) (see [15]).

B. Deep LSTM-RNN
An LSTM-RNN is a neural network designed for processing sequential data with inputs {x i } N i=0 and outputs {y i } N i=0 , and defined as The activation function φ is given by the following relations: where σ is the logistic sigmoid function and W and b terms represent weight matrices and bias vectors to be optimized, respectively. The deep LSTM-RNN can be constructed by stacking multiple of these layers [1], [2].
Since contraction analysis for discrete-time systems leads to similar results [5], [11], we define the inputs x i as discretized states {x i = x(t i )} N t=0 , and the outputs y i as non-zero components of the unique Cholesky decomposition of the optimal contraction metric, as will be discussed in Sec. III.

III. NEURAL CONTRACTION METRIC (NCM)
This section presents an algorithm to obtain an NCM depicted in Fig. 1.

A. Convex Optimization-based Sampling of Contraction Metrics (CV-STEM)
We derive one approach to sample contraction metrics by using the ConVex optimization-based Steady-state Tracking Error Minimization (CV-STEM) method [10], [11], which could handle the control design of Itô stochastic nonlinear systems. In this section, we propose its simpler formulation for nonlinear systems with bounded disturbances in order to be of practical use in engineering applications.
By Theorem 1, the problem to minimize an upper bound of the steady-state Euclidean distance between the trajectory x 1 (t) of the unperturbed system and x 2 (t) of the perturbed system (1) ((4) as t → ∞) can be formulated as follows: (2) and (3)  (5) where W (x,t) = M(x,t) −1 is used as a decision variable instead of M(x,t). We assume that the contraction rate α and disturbance bound d are given in (5) (see Remark 1 on how to select α). We need the following lemma to convexify this nonlinear optimization problem. Lemma 1: The inequalities (2) and (3) are equivalent tȯ respectively, where χ = ω/ω,W = νW , and ν = 1/ω. Proof: Since ν = 1/ω > 0 and W (x,t) 0, multiplying (2) by ν and then by W (x,t) from both sides preserves matrix definiteness and the resultant inequalities are equivalent to the original ones [27, pp. 114]. These operations yield (6). Next, since M(x,t) 0, there exists a unique µ(x,t) 0 s.t. M = µ 2 . Multiplying (3) by µ −1 from both sides gives ωI W (x,t) ωI as we have ω, ω > 0 and (µ −1 ) 2 = W . We get (7) by multiplying this inequality by ν = 1/ω. We are now ready to state and prove our main result on the convex optimization-based sampling.
Remark 1: Since (6) and (7) are independent of ν = 1/ω, the choice of ν does not affect the optimal value of the minimization problem in Theorem 2. In practice, as we have sup x,t M(x,t) ≤ 1/ω = ν by (3), it can be used as a penalty to optimally adjust the induced 2-norm of estimation and control gains when the problem explicitly depends on ν (see Sec. IV for details). Also, although α is fixed in Theorem 2, it can be found by a line search as will be demonstrated in Sec. V.
Remark 2: The problem (8) can be solved as a finitedimensional problem by using backward difference ap- ∀i with ∆t ∆t 2 > 0, and by discretizing it along a pre-computed system trajectory {x(t i )} N i=0 .

B. Deep LSTM-RNN Training
Instead of directly using sequential data of optimal contraction metrics {M(x(t i ),t i )} N i=0 for neural network training, the positive definiteness of M(x,t) is utilized to reduce the dimension of the target output {y i } N i=0 defined in Sec. II. . We denote these nonzero entries as θ (x,t) ∈ R 1 2 n(n+1) . As a result, the dimension of the target data θ (x,t) is reduced by n(n − 1)/2 without losing any information on M(x,t).
The pseudocode to obtain an NCM depicted in Fig. 1 is presented in Algorithm 1. The deep LSTM-RNN in Sec. II is trained with the sequential state data {x(t i )} N i=0 and the target data {θ (x(t i ),t i )} N i=0 using Stochastic Gradient Descent (SGD). We note these pairs will be sampled for multiple trajectories to increase sample size and to avoid overfitting.
using SGD Compute the test error for data in S test if test error < ε then break

IV. NCM-BASED OPTIMAL ESTIMATION AND CONTROL
This section delineates how to construct an NCM offline and utilize it online for state estimation and feedback control.

A. Problem Statement
We apply an NCM to the state estimation problem for the following nonlinear system with bounded disturbances: , and C(x,t) = (∂ h/∂ x). We design an estimator aṡ W +WA(x,t) + A(x,t) T W − 2C(x,t) T C(x,t) −2αW (11) where α > 0. The virtual system of (9) and (10) is given aṡ where d e (q,t) is defined as d e (x,t) = B(x,t)d 1 (t) and d e (x,t) = M(x,t)C(x,t) T G(x,t)d 2 (t). Note that (12) has q = x and q =x as its particular solutions. The differential dynamics of (12) with d e = 0 is given as

B. Nonlinear Stability Analysis
We have the following lemma for deriving a condition to guarantee the local contraction of (12) in Theorem 3.
Proof: See Lemma 2 of [29] or Theorem 1 of [9]. The following theorem along with this lemma guarantees the exponential stability of the estimator (10).

C. Convex Optimization-based Sampling (CV-STEM)
We have the following proposition to sample optimal contraction metrics for the NCM-based state estimation.
We have an analogous result for state feedback control.
Corollary 1: Consider the following system and a state feedback controller u(t) with the bounded disturbance d(t): , assuming that f (x,t) = 0 at x = 0 [24], [25]. Suppose there exist positive constants ω, ω, and b 2 s.t. ωI W (x,t) ωI and B 2 (x,t) ≤ b 2 , ∀x,t. Then M(x,t) that minimizes an upper bound of lim t→∞ x 0 δ q can be found by the following convex optimization problem: where χ = ω/ω, ν = 1/ω,W = νW , and λ > 0 is a userdefined constant. The arguments of A(x,t), B 1 (x,t), andW (x,t) are omitted for notational simplicity. Proof: The system with q = x, 0 as its particular solutions is given byq where d c (x,t) = B 2 (x,t)d(t) and d c (0,t) = 0. Since we have d c (q,t) ≤ b 2 d and the differential dynamics is when d c = 0, we get lim t→∞ x 0 δ q ≤ b 2 dχ/α by the same proof as for Theorem 3 with (13) replaced by (20), (14) by (18) then follows as in the proof of Proposition 1, where λ ≥ 0 is for penalizing excessively large control inputs through ν ≥ sup x,t M(x,t) (see Remark 1).

D. NCM Construction and Interpretation
Algorithm 1 along with Proposition 1 and Corollary 1 returns NCMs to computex(t) of (10) and u(t) of (17) for state estimation and control in real-time. They also provide the bounded error tube (see Theorem 1, [15], [16]) for robust motion planning problems as will be seen in Sec. V.
The similarity of Corollary 1 to Proposition 1 stems from the estimation and control duality due to the differential nature of contraction analysis as is evident from (13) and (20). Analogously to the discussion of the Kalman filter and LQR duality in LTV systems, this leads to two different interpretations on the weight of ν (d 2cḡ /γ in (16) and λ in (19)). As discussed in Remark 1, one way is to see it as a penalty on the induced 2-norm of feedback gains. Since d 2 = 0 in (16) means no noise acts on y(t), it can also be viewed as an indicator of how much we trust the measurement y(t): the larger the weight of ν, the less confident we are in y(t). These agree with our intuition as smaller feedback gains are suitable for systems with larger measurement uncertainty.

V. SIMULATION
The NCM framework is demonstrated using Lorentz oscillator state estimation and spacecraft motion planning and control problems. CVXPY [13] with the MOSEK solver [14] is used to solve convex optimization problems. A Python implementation of Algorithm 1 is available at https://github.com/astrohiro/ncm.

A. State Estimation of Lorenz Oscillator
We first consider state estimation of the Lorentz oscillator with bounded disturbances described asẋ = f (x) + d 1 (t) and , and sup t d 2 (t) = 1. We use dt = 0.1 for integration, with one measurement y per dt.
1) Sampling of Optimal Contraction Metrics: Using Proposition 1, we sample the optimal contraction metric along 100 trajectories with uniformly distributed initial conditions (−10 ≤ x i ≤ 10, i = 1, 2, 3). Figure 2 plots J * CVe in (16) for 100 different trajectories and the optimal α is found to be α = 3.4970. The optimal estimator parameters averaged over 100 trajectories for α = 3.4970 are summarized in Table I. 2) Deep LSTM-RNN Training: A deep LSTM-RNN is trained using Algorithm 1 and Proposition 1 with the sequential data {{(x s (t i ), θ s (x(t i ))} N i=0 } S s=1 sampled over the 100 different trajectories (S = 100). Note that θ s (x(t i )) are standardized and normalized to make the SGD-based learning process stable. Figure 3 shows the test loss of the LSTM-RNN models with different number of layers and hidden units. We can see that the models with more than 2 layers overfit and those with less than 32 hidden units underfit to the training samples. Thus, the number of layers and hidden units are selected as 2 and 64, respectively. The resultant MSE of the trained LSTM-RNN is shown in Table I.
3) State Estimation with an NCM: The estimation problem is solved using the NCM, sampling-based CV-STEM [10], [11], and Extended Kalman Filter (EKF) with sup t d 1 (t) = 20    while the EKF has a larger error compared to the other two. As expected from the small MSE of Table I, the estimation error of the NCM estimator shows a trend similar to that of the sampling-based CV-STEM estimator without losing its estimation performance.

B. Spacecraft Optimal Motion Planning
We consider an optimal motion planning problem of the planar spacecraft dynamical system, given asẋ = Ax + B(x)u + d(t), where u ∈ R 8 , sup t d(t) = 0.15, and x = [p x , p y , φ ,ṗ x ,ṗ y ,φ ] T with p x , p y , and φ being the horizontal coordinate, vertical coordinate, and yaw angle of the spacecraft, respectively. The constant matrix A and the statedependent actuation matrix B(x) are defined in [30]. All the parameters of the spacecraft are normalized to 1.
2) Motion Planning with an NCM: Given an NCM, we can solve a robust motion planning problem, where the state constraint is now described by the bounded error tube (see Theorem 1) with radius R tube = d( √ χ/α) = 0.4488. Figure 5 parameters ν χ J *   shows the spacecraft motion (p x , p y ) on a planar field, computed using the NCM, sampling-based CV-STEM [10], [11], and baseline LQR control with Q = 2.4I and R = I which does not account for the disturbance. As summarized in Table III, input constraints 0 ≤ u i (t) ≤ 1, ∀i, ∀t are satisfied and the three controllers use a similar amount of control effort. All the controllers except the LQR keep their trajectories within the tube avoiding collision with the circular obstacles, even under the presence of disturbances as depicted in Fig. 5. Also, the NCM controller predicts the computationally-expensive CV-STEM controller with the small MSE as given in the last column of Table II. VI. CONCLUSION In this paper, we present a Neural Contraction Metric (NCM), a deep learning-based global approximation of an optimal contraction metric for online nonlinear estimation and control. The novelty of the NCM approach lies in that: 1) data points of the optimal contraction metric are sampled offline by solving a convex optimization problem, which minimizes an upper bound of the steady-state Euclidean distance between the perturbed and unperturbed trajectories without assuming any hypothesis function space, and 2) the deep LSTM-RNN is constructed to model the sampled metrics with arbitrary accuracy. The superiority of its performance is validated through numerical simulations on state estimation and optimal motion planning problems.