IEEE CONTROL SYSTEMS LETTERS, VOL. 5, NO. 1, JANUARY 2021 211
Neural Contraction Metrics for Robust Estimation
and Control: A Convex Optimization Approach
Hiroyasu Tsukamoto , Graduate Student Member, IEEE, and Soon-Jo Chung , Senior Member, IEEE
Abstract—This letter presents a new deep learning-
based framework for robust nonlinear estimation and con-
trol 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 non-
linear 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 per-
turbed and unperturbed system trajectories. We demon-
strate how to exploit NCMs to design an online optimal esti-
mator 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.
Index Terms—Machine learning, observers for nonlinear
systems, optimal control.
I. INTRODUCTION
P
ROVABLY stable and optimal state estimation and con-
trol 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 con-
traction 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 prov-
ably 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
Manuscript received March 17, 2020; revised May 15, 2020; accepted
June 4, 2020. Date of publication June 11, 2020; date of current
version June 24, 2020. This work was supported in part by the Jet
Propulsion Laboratory, California Institute of Technology and in part by
the Raytheon Company. Recommended by Senior Editor G. Cherubini.
(Corresponding author: Hiroyasu Tsukamoto.)
The authors are with the Graduate Aerospace Laboratories,
California Institute of Technology, Pasadena, CA 91125 USA (e-mail:
htsukamoto@caltech.edu; sjchung@caltech.edu).
Data is available on-line at https://github.com/astrohiro/ncm
Digital Object Identifier 10.1109/LCSYS.2020.3001646
improved memory structure proposed to circumvent gradient
vanishing [3] and is a universal approximator of contin-
uous 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 func-
tion space. These sampled metrics, the existence of which is
a necessary and sufficient condition for exponential conver-
gence [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 dynam-
ics [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 per-
turbed and unperturbed system trajectories. In this letter, we
present a convex optimization problem equivalent to this
problem, thereby exploiting the differential nature of con-
traction 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 effi-
cient 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 feed-
back estimation and control gain or a bounded error tube for
robust motion planning [15], [16].
We illustrate how to design an optimal NCM-based esti-
mator and controller for nonlinear systems with bounded
disturbances, utilizing the estimation and control duality in dif-
ferential 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
This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/
212 IEEE CONTROL SYSTEMS LETTERS, VOL. 5, NO. 1, JANUARY 2021
Fig. 1. Illustration of the NCM: M(x, t) denotes the optimal contrac-
tion metric; x(t) and x
d
(t) denote perturbed and unperturbed system
trajectories; h
i
and c
i
denote the hidden states of the deep LSTM-RNN,
respectively.
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 candi-
date 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 con-
ditions. 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 nonlin-
ear 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. P
RELIMINARIES
We use the notation x and A for the Euclidean and
induced 2-norm, and A 0, A 0, A ≺ 0, and A 0for
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:
˙x(t) = f (x(t), t) + d(t) (1)
where t ∈ R
≥0
, x : R
≥0
→ R
n
, f : R
n
× R
≥0
→ R
n
, and
d : R
≥0
→ R
n
with d = sup
t≥0
d(t) < +∞.
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
M(x, t) = (x, t)
T
(x, t) 0, α>0, and 0 <ω, ω<∞ s.t.
˙
M(x, t) + sym
M(x, t)
∂f(x, t)
∂x
−2αM(x, t), ∀x, t (2)
ω
−1
I M(x, t) ω
−1
I, ∀x, t. (3)
Then the smallest path integral is exponentially bounded,
thereby yielding a bounded tube of states:
x
2
x
1
δx−R(0)
√
ωe
−αt
≤
d
α
ω
ω
≤
dω
αω
(4)
where R(t) =
x
2
x
1
(x, t)δx(t), ∀t.
Proof: Differentiating R(t) yields
˙
R(t) + αR(t)
≤(x(t), t)d(t) (see [5]). Since we have (x(t), t)d(t)
≤
d/
√
ω, applying the comparison lemma [26] results in
R(t) ≤ R(0)e
−αt
+ d/(α
√
ω). 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 process-
ing sequential data with inputs {x
i
}
N
i=0
and outputs {y
i
}
N
i=0
,
and defined as y
i
= W
hy
h
i
+ b
y
where h
i
= φ(x
i
, h
i−1
, c
i−1
).
The activation function φ is given by the following relations:
h
i
= o
i
tanh(c
i
), o
i
= σ(W
xo
x
i
+ W
ho
h
i−1
+ W
co
c
i
+ b
o
), c
i
=
f
i
c
i−1
+ι
i
tanh(W
xc
x
i
+W
hc
h
i−1
+b
c
), f
i
= σ(W
xf
x
i
+W
hf
h
i−1
+
W
cf
c
i−1
+b
f
), and ι
i
= σ(W
xi
x
i
+W
hi
h
i−1
+W
ci
c
i−1
+b
i
), 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 compo-
nents of the unique Cholesky decomposition of the optimal
contraction metric, as will be discussed in Section III.
III. N
EURAL 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) (i.e., (4) as t →∞) can be formulated as follows:
J
∗
NL
= min
ω>0,ω>0,W0
dω
αω
s.t. (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 aregivenin(5)(seeRemark1onhow
to select α). We need the following lemma to convexify this
nonlinear optimization problem.
TSUKAMOTO AND CHUNG: NEURAL CONTRACTION METRICS FOR ROBUST ESTIMATION AND CONTROL: CONVEX OPTIMIZATION APPROACH 213
Lemma 1: The inequalities (2) and (3) are equivalent to
˙
˜
W(x, t) − sym
∂f (x, t)
∂x
˜
W(x, t)
2α
˜
W(x, t), ∀x, t (6)
I
˜
W(x, t) χI, ∀x, t (7)
respectively, where χ =
ω/ω,
˜
W = νW, and ν = 1/ω.
Proof: Since ν = 1/ω
> 0 and W(x, t) 0, multiply-
ing (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, p. 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.
Theorem 2: Consider the convex optimization problem:
J
∗
CV
= min
χ∈R,
˜
W0
dχ
α
s.t. (6) and (7) (8)
where χ and
˜
W are defined in Lemma 1, and α>0 and
d = sup
t≥0
d(t) are assumed to be given. Then J
∗
NL
= J
∗
CV
.
Proof: By definition, we have
dω/(αω) = dχ/α. Since (2)
and (3) are equivalent to (6) and (7) by Lemma 1, rewriting
the objective in the original problem (5) using this equality
completes the proof.
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 con-
trol gains when the problem explicitly depends on ν (see
Section IV for details). Also, although α is fixed in Theorem 2,
it can be found by a line search as will be demonstrated in
Section V.
Remark 2: The problem (8) can be solved as a finite-
dimensional problem by using backward difference approx-
imation,
˙
˜
W(x(t
i
), t
i
) (
˜
W(x(t
i
), t
i
) −
˜
W(x(t
i−1
), t
i−1
))/t
i
,
where t = t
i
− t
i−1
, ∀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 con-
traction 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 Section II.
Lemma 2: AmatrixA 0 has a unique Cholesky decom-
position, i.e., there exists a unique upper triangular matrix
U ∈ R
n×n
with strictly positive diagonal entries s.t. A = U
T
U.
Proof: See [28, p. 441].
As a result of Lemma 2, we can select (x, t) defined in
Theorem 1 as the unique Cholesky decomposition of M(x, t)
and train the deep LSTM-RNN using only the non-zero entries
of the unique upper triangular matrices {(x(t
i
), t
i
)}
N
i=0
.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 Section II
is trained with the sequential state data {x(t
i
)}
N
i=0
and the target
Algorithm 1: NCM Algorithm
Inputs : Initial and terminal states {x
s
(t
0
), x
s
(t
N
)}
S
s=1
Outputs: NCM and steady-state bound J
∗
CV
in (8)
A. Sampling of Optimal Contraction Metrics
for s ← 1 to S do
Generate a trajectory {x
s
(t
i
)}
N
i=0
using x
s
(t
0
) (could use x
s
(t
N
) for
motion planning problems)
for α ∈ A
linesearch
do
Find J
∗
CV
(α, x
s
) and {θ(x
s
(t
i
), t
i
)}
N
i=0
by Th. 2
Find α
∗
(x
s
) = arg min
α∈A
linesearch
J
∗
CV
(α, x
s
)
Save {θ(x
s
(t
i
), t
i
)}
N
i=0
for α = α
∗
(x
s
)
Obtain J
∗
CV
= max
s
J
∗
CV
(α
∗
(x
s
), x
s
)
B. Deep LSTM-RNN Training
Split data into a train set S
train
and test set S
test
for epoch ← 1 to N
epochs
do
for s ∈ S
train
do
Train the deep LSTM-RNN with {x
s
(t
i
)}
N
i=0
, {θ(x
s
(t
i
), t
i
)}
N
i=0
using SGD
Compute the test error for data in S
test
if test error < then
break
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.
IV. NCM-B
ASED 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:
˙x = f (x, t) + B(x, t)d
1
(t), y(t) = h(x, t) + G(x, t)d
2
(t) (9)
where d
1
: R
≥0
→ R
k
1
, B : R
n
× R
≥0
→ R
n×k
1
,
y : R
≥0
→ R
m
, d
2
: R
≥0
→ R
k
2
, h : R
n
× R
≥0
→ R
m
,
and G : R
n
× R
≥0
→ R
m×k
2
with d
1
= sup
t
d
1
(t) < +∞
and d
2
= sup
t
d
2
(t) < +∞.LetW = M(ˆx, t)
−1
0,
A(x, t) = (∂f /∂x), and C(x, t) = (∂h/∂x).Letˆx : R
≥0
→ R
n
.
We design an estimator as
˙
ˆx = f (ˆx, t) + M(ˆx, t)C(ˆx, t)
T
(y − h(ˆx, t)) (10)
˙
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 as
˙q = f (q, t) + M(ˆx, t)C(ˆx, t)
T
(h(x, t) − h(q, t)) + d
e
(q, t) (12)
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
δ ˙q = (A(q, t) − M(ˆx, t)C(ˆx, t)
T
C(q, t))δq. (13)
B. Nonlinear Stability Analysis
We have the following lemma for deriving a condition to
guarantee the local contraction of (12) in Theorem 3.
Lemma 3: If (11) holds for t ≥ 0, there exists r(t)>0 s.t.
2γ W +
˙
W + sym (WA(q, t)) − sym(C(ˆx, t)
T
C(q, t)) 0 (14)
214 IEEE CONTROL SYSTEMS LETTERS, VOL. 5, NO. 1, JANUARY 2021
for all q(t) with q(t) −ˆx(t)≤r(t), where 0 <γ <α.
Proof: See [29, Lemma 2] or [9, Theorem 1].
The following theorem along with this lemma guarantees
the exponential stability of the estimator (10).
Theorem 3: Suppose that there exist positive constants ω
,
ω, b, ¯c, ¯g, and ρ s.t. ωI W(ˆx, t) ωI, B(x, t)≤b,
C(ˆx, t)≤¯c, G(x, t)≤¯g, and r(t) ≥ ρ, ∀ˆx, x, t, where r(t)
is defined in Lemma 3. If (11) holds and R
e
(0)+D
e
/γ ≤
√
ωρ,
where R
e
(t) =
x
ˆx
(ˆx, t)δq(t) with W =
T
and D
e
=
d
1
b
√
ω + d
2
¯c¯g/
√
ω, then the distance between the trajectory
of (9) and (10) is exponentially bounded as follows:
x
ˆx
δq≤
R
e
(0)
√
ω
e
−γ t
+
d
1
b
γ
χ +
d
2
¯c¯g
γ
ν (15)
where χ =
ω/ω, ν = 1/ω, and 0 <γ <α.
Proof: Using (13), we have d(δq
2
)/dt = δq
T
(
˙
W +
sym(WA(q, t)) − sym(C(ˆx, t)
T
C(q, t)))δq when d
e
= 0. This
along with (14) gives
˙
R
e
(t) ≤−γ R
e
(t) in the region
where Lemma 3 holds. Thus, using the bound (ˆx(t), t)d
e
(q, t)≤D
e
,wehave
√
ω
x
ˆx
δq≤R
e
(0)e
−γ t
+ D
e
/γ
by the same proof as for Theorem 1. Rewriting this with
χ, ν, and 1 ≤
√
χ ≤ χ yields (15). This also implies that
√
ωx −ˆx≤R
e
(0) + D
e
/γ, ∀t. Hence, the sufficient condi-
tion for q −ˆx in Lemma 3 reduces to the one required in
this theorem.
C. Convex Optimization-Based Sampling (CV-STEM)
We have the following proposition to sample optimal con-
traction metrics for the NCM-based state estimation.
Proposition 1: M(ˆx, t) that minimizes an upper bound of
lim
t→∞
x
ˆx
δq is found by the convex optimization problem:
J
∗
CVe
= min
ν>0,χ ∈R,
˜
W0
d
1
b
γ
χ +
d
2
¯c¯g
γ
ν
s.t.
˙
˜
W +
˜
WA + A
T
˜
W − 2νC
T
C −2α
˜
W and I
˜
W χI (16)
where χ =
ω/ω, ν = 1/ω,
˜
W = νW, and 0 <γ <α.
The arguments of A(ˆx, t), C(ˆx, t), and
˜
W(ˆx, t) are omitted for
notational simplicity.
Proof: Multiplying (11) and ω
I W(ˆx, t) ωI, ∀ˆx, t by ν
yields the constraints of (16). Then applying Theorem 2 with
the objective function given in (15) of Theorem 3 as t →∞
yields (16).
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):
˙x = f (x, t) + B
1
(x, t)u + B
2
(x, t)d(t) (17)
˙
W − A(x, t)W − WA(x, t)
T
+ 2B
1
(x, t)B
1
(x, t)
T
2αW (18)
where u =−B
1
(x, t)
T
M(x, t)x, B
1
: R
n
× R
≥0
→ R
n×m
,
B
2
: R
n
× R
≥0
→ R
n×k
, W = M
−1
0, α>0, and A is a
matrix defined as A(x, t)x = f(x, 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:
J
∗
CVc
= min
ν>0,χ ∈R,
˜
W0
b
2
d
α
χ + λν
s.t. −
˙
˜
W + A
˜
W +
˜
WA
T
− 2νB
1
B
T
1
−2α
˜
W and I
˜
W χ I (19)
where χ =
ω/ω, ν = 1/ω,
˜
W = νW, and λ>0isa
user-defined constant. The arguments of A(x, t), B
1
(x, t), and
˜
W(x, t) are omitted for notational simplicity.
Proof: The system with q = x, 0 as its particular solutions
is given by ˙q = (A(x, t) −B
1
(x, t)B
1
(x, t)
T
M(x, t))q+d
c
(q, t),
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
δ ˙q = (A(x, t) − B
1
(x, t)B
1
(x, t)
T
M(x, t))δq (20)
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), and R
e
(t) by R
c
(t) =
x
0
(x, t)δq(t), where
M =
T
. Equation (19) 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 compute ˆx(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 Section 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 interpreta-
tions on the weight of ν (i.e.,
d
2
¯c¯g/γ 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
= 0in(16)
means no noise acts on y(t), it can also be viewed as an indi-
cator 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
measurements with larger uncertainty.
V. S
IMULATION
The NCM framework is demonstrated using Lorentz oscilla-
tor 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
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 ˙x = f (x) + d
1
(t) and
y = Cx+d
2
(t), where f (x) = [σ(x
2
−x
1
), x
1
(ρ−x
3
)−x
2
, x
1
x
2
−
βx
3
]
T
, x = [x
1
, x
2
, x
3
]
T
, σ = 10, β = 8/3, ρ = 28, C =
[100],sup
t
d
1
(t)=
√
3, 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 sequen-
tial data {{(x
s
(t
i
), θ
s
(x(t
i
))}
N
i=0
}
S
s=1
sampled over the 100 differ-
ent trajectories (S = 100). Note that θ
s
(x(t
i
)) are standardized
TSUKAMOTO AND CHUNG: NEURAL CONTRACTION METRICS FOR ROBUST ESTIMATION AND CONTROL: CONVEX OPTIMIZATION APPROACH 215
Fig. 2. Upper bound of steady-state errors as a function of α: each
curve with a different color corresponds to a trajectory with a different
initial condition.
TABLE I
NCM E
STIMATOR PARAMETERS FOR α = 3.4970 WITH J
∗
CVe
AND MSE
A
VERAGED OVER 100 TRAJECTORIES
Fig. 3. LSTM-RNN test loss with 2 layers (left) and 32 hidden
units (right).
and normalized to make the SGD-based learning process sta-
ble. 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
√
3,
and sup
t
d
2
(t)=20. We use x(0) = [−1.0, 2.0, 3.0]
T
and
ˆx(0) = [150.1, −1.5, −6]
T
for the actual and estimated initial
conditions, respectively. The EKF weight matrices are selected
as R = 20I and Q = 10I.
Figure 4 shows the smoothed estimation error x(t) −ˆx(t)
using a 15-point moving average filter. The errors of the
NCM and CV-STEM estimators are below the optimal upper
bound while the EKF has a larger error compared to the
other two. As expected from the small MSE of Table I,
Figure 4 implies that the NCM estimator is able to approxi-
mate 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 ˙x = Ax +
B(x)u + d(t), where u ∈ R
8
,sup
t
d(t)=0.15, and
Fig. 4. Lorentz oscillator state estimation error smoothed using a
15-point moving average filter.
TABLE II
NCM C
ONTROL PARAMETERS FOR α = 0.58 WITH J
∗
CVc
AND MSE
A
VERAGED OVER 100 TRAJECTORIES
TABLE III
C
ONTROL PERFORMANCES FOR SPACECRAFT
MOTION PLANNING PROBLEM
x = [p
x
, p
y
,φ,˙p
x
, ˙p
y
,
˙
φ]
T
with p
x
, p
y
, and φ being the hor-
izontal coordinate, vertical coordinate, and yaw angle of the
spacecraft, respectively. The constant matrix A and the state-
dependent actuation matrix B(x) are defined in [30]. All the
parameters of the spacecraft are normalized to 1.
1) Problem Formulation: In the planar field, we have 6 cir-
cular obstacles with radius 3 located at (p
x
, p
y
) = (0, 11),
(5, 3), (8, 11), (13, 3), (16, 11), and (21, 3). The goal of the
motion planning problem is to find an optimal trajectory that
avoids these obstacles and minimize
50
0
u(t)
2
dt subject to
input constraints 0 ≤ u
i
(t) ≤ 1, ∀i, ∀t and the dynamics con-
straints. The initial and terminal condition are selected as
x(0) = [0, 0,π/12, 0, 0, 0]
T
and x(t
N
) = [20, 18, 0, 0, 0, 0]
T
.
Following the same procedure described in the state estima-
tion problem, the optimal control parameters and the MSE of
the LSTM-RNN trained using Algorithm 1 with Corollary 1
are determined as shown in Table II.
2) Motion Planning With an NCM: GivenanNCM,wecan
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
shows the spacecraft motion (p
x
, p
y
) on a planar field, com-
puted 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
control effort of these controllers are adjusted to be below 50
for the sake of fair comparison. All the controllers except the
LQR keep their trajectories within the tube avoiding collision
with the circular obstacles, even under the presence of distur-
bances as depicted in Fig. 5. Also, the NCM controller predicts
