IEEE CONTROL SYSTEMS LETTERS (L-CSS), PREPRINT VERSION, ACCEPTED JUNE, 2020.
1
Neural Contraction Metrics for Robust Estimation
and Control: A Convex Optimization Approach
Hiroyasu Tsukamoto,
Member, IEEE
, and Soon-Jo Chung,
Senior Member, IEEE
Abstract
—This paper presents a new deep learning-based
framework for robust nonlinear estimation and control using
the concept of a N
eural C
ontraction M
etric (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 ex-
ploit 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.
Index Terms
—Machine learning, Observers for nonlinear sys-
tems, Optimal control.
I. I
NTRODUCTION
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 contrac-
tion 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 as-
suming 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
This work was in part funded by the Jet Propulsion Laboratory, California
Institute of Technology and the Raytheon Company.
The authors are with the Graduate Aerospace Laboratories, California
Institute of Technology, Pasadena, CA, USA. E-mail:
{
htsukamoto,
sjchung
}
@caltech.edu
), Code: https://github.com/astrohiro/ncm.
Fig. 1: Illustration of the NCM:
M
(
x
,
t
)
denotes the optimal
contraction metric;
x
(
t
)
and
x
d
(
t
)
denote perturbed and unper-
turbed system trajectories;
h
i
and
c
i
denote the hidden states
of the deep LSTM-RNN, respectively.
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 up-
per bound of the steady-state Euclidean distance between
perturbed and unperturbed system trajectories. In this pa-
per, 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
ˆ
o 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 es-
timator 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
arXiv:2006.04361v1 [eess.SY] 8 Jun 2020
2
IEEE CONTROL SYSTEMS LETTERS (L-CSS), PREPRINT VERSION, ACCEPTED JUNE, 2020.
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 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
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. P
RELIMINARIES
We use the notation
‖
x
‖
and
‖
A
‖
for the Euclidean and in-
duced 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:
̇
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
)
6
=
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 rela-
tions:
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
x f
x
i
+
W
h f
h
i
−
1
+
W
c f
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 Sec. III.
III. N
EURAL
C
ONTRACTION
M
ETRIC
(NCM)
This section presents an algorithm to obtain an NCM
depicted in Fig. 1.
A. Convex Optimization-based Sampling of Contraction Met-
rics (CV-STEM)
We derive one approach to sample contraction metrics by
using the C
onV
ex optimization-based S
teady-state T
racking
E
rror M
inimization (CV-STEM) method [10], [11], which
could handle the control design of It
ˆ
o 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:
J
∗
NL
=
min
ω
>
0
,
ω
>
0
,
W
0
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
are given in (5) (see Remark 1 on how
H. TSUKAMOTO
et al.
: NEURAL CONTRACTION METRICS FOR ROBUST ESTIMATION AND CONTROL: A CONVEX OPTIMIZATION APPROACH
3
to select
α
). We need the following lemma to convexify this
nonlinear optimization problem.
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, 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.
Theorem 2:
Consider the convex optimization problem:
J
∗
CV
=
min
χ
∈
R
,
̃
W
0
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 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 finite-
dimensional problem by using backward difference ap-
proximation,
̇
̃
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 Sec. II.
Lemma 2:
A matrix
A
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, pp. 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 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.
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
IV. NCM-B
ASED
O
PTIMAL
E
STIMATION AND
C
ONTROL
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
)
‖
<
+
∞
. Let
W
=
M
(
ˆ
x
,
t
)
−
1
0,
A
(
x
,
t
) = (
∂
f
/
∂
x
)
, and
C
(
x
,
t
) =
(
∂
h
/
∂
x
)
. 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
−
2
C
(
ˆ
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)
4
IEEE CONTROL SYSTEMS LETTERS (L-CSS), PREPRINT VERSION, ACCEPTED JUNE, 2020.
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)
for all
q
(
t
)
with
‖
q
(
t
)
−
ˆ
x
(
t
)
‖≤
r
(
t
)
, where 0
<
γ
<
α
.
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).
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
, we have
√
ω
∫
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 condition 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
,
̃
W
0
d
1
b
γ
χ
+
d
2
̄
c
̄
g
γ
ν
(16)
s.t.
̇
̃
W
+
̃
WA
+
A
T
̃
W
−
2
ν
C
T
C
−
2
α
̃
W
and
I
̃
W
χ
I
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 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
+
2
B
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
,
̃
W
0
b
2
d
α
χ
+
λ ν
(19)
s.t.
−
̇
̃
W
+
A
̃
W
+
̃
WA
T
−
2
ν
B
1
B
T
1
−
2
α
̃
W
and
I
̃
W
χ
I
where
χ
=
ω
/
ω
,
ν
=
1
/
ω
,
̃
W
=
ν
W
, and
λ
>
0 is a 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
Θ
.
(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 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
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
=
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. 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
of Algorithm 1 is available at https://github.com/astrohiro/ncm.
H. TSUKAMOTO
et al.
: NEURAL CONTRACTION METRICS FOR ROBUST ESTIMATION AND CONTROL: A CONVEX OPTIMIZATION APPROACH
5
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.
parameters
ν
χ
J
∗
CVe
MSE
values
1
.
3375
×
10
2
9
.
2977
42
.
852
1
.
0190
×
10
−
3
TABLE I: NCM estimator parameters for
α
=
3
.
4970:
J
∗
CVe
and the MSE are averaged over 100 trajectories.
A. State Estimation of Lorenz Oscillator
We first consider state estimation of the Lorentz oscilla-
tor 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
= [
1 0 0
]
, 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 Propo-
sition 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 se-
quential 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 stan-
dardized 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
√
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
=
20
I
and
Q
=
10
I
.
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
Fig. 3: LSTM-RNN test loss with 2 layers (left) and 32 hidden
units (right).
Fig. 4: Lorentz oscillator state estimation error smoothed using
a 15-point moving average filter.
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 ̇
x
=
Ax
+
B
(
x
)
u
+
d
(
t
)
, where
u
∈
R
8
, sup
t
‖
d
(
t
)
‖
=
0
.
15, and
x
=
[
p
x
,
p
y
,
φ
,
̇
p
x
,
̇
p
y
,
̇
φ
]
T
with
p
x
,
p
y
, and
φ
being the horizon-
tal 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
circular 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
constraints. 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 estimation
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:
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
6
IEEE CONTROL SYSTEMS LETTERS (L-CSS), PREPRINT VERSION, ACCEPTED JUNE, 2020.
parameters
ν
χ
J
∗
CVc
MSE
values
6
.
2090
×
10
2
3
.
0116
8
.
9524
3
.
4675
×
10
−
4
TABLE II: NCM control parameters for
α
=
0
.
58:
J
∗
CVc
and
the MSE are averaged over 100 trajectories.
Fig. 5: Spacecraft motion
(
p
x
,
p
y
)
on a planar field.
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
.
4
I
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. C
ONCLUSION
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.
∫
50
0
‖
u
(
t
)
‖
2
dt
min
i
(
inf
t
u
i
(
t
))
max
i
(
sup
t
u
i
(
t
))
NCM
4
.
8959
×
10
−
1
.
2728
×
10
−
8
1
.
0000
CV-STEM
4
.
8019
×
10
−
4
.
1575
×
10
−
8
1
.
0000
LQR
4
.
9260
×
10
0
.
0000
1
.
0000
TABLE III: Control performances for the spacecraft motion
planning problem.
R
EFERENCES
[1] S. Hochreiter and J. Schmidhuber, “Long short-term memory,”
Neural
Computation
, vol. 9, no. 8, pp. 1735–1780, 1997.
[2] A. Graves, A. Mohamed, and G. Hinton, “Speech recognition with deep
recurrent neural networks,” in
IEEE Int. Conf. Acoust. Speech Signal
Process.
, May 2013, pp. 6645–6649.
[3] Y. Bengio, P. Simard, and P. Frasconi, “Learning long-term dependencies
with gradient descent is difficult,”
IEEE Trans. Neural Networks
, vol. 5,
no. 2, pp. 157–166, Mar. 1994.
[4] K. Funahashi and Y. Nakamura, “Approximation of dynamical systems
by continuous time recurrent neural networks,”
Neural Networks
, vol. 6,
no. 6, pp. 801 – 806, 1993.
[5] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for nonlinear
systems,”
Automatica
, vol. 34, no. 6, pp. 683 – 696, 1998.
[6] G. Shi
et al.
, “Neural lander: Stable drone landing control using learned
dynamics,” in
IEEE Int. Conf. Robot. Automat.
, May 2019.
[7] S. M. Richards, F. Berkenkamp, and A. Krause, “The Lyapunov neural
network: Adaptive stability certification for safe learning of dynamical
systems,” in
Conf. Robot Learn.
, vol. 87, Oct. 2018, pp. 466–476.
[8] Y.-C. Chang, N. Roohi, and S. Gao, “Neural Lyapunov control,” in
Adv.
Neural Inf. Process. Syst.
, 2019, pp. 3245–3254.
[9] A. P. Dani, S.-J. Chung, and S. Hutchinson, “Observer design for
stochastic nonlinear systems via contraction-based incremental stability,”
IEEE Trans. Autom. Control
, vol. 60, no. 3, pp. 700–714, Mar. 2015.
[10] H. Tsukamoto and S.-J. Chung, “Convex optimization-based controller
design for stochastic nonlinear systems using contraction analysis,” in
58th IEEE Conf. Decis. Control
, Dec. 2019, pp. 8196–8203.
[11] ——, “Robust controller design for stochastic nonlinear systems via
convex optimization,”
Submitted to IEEE Trans. Autom. Control
, 2019.
[12] S. Boyd and L. Vandenberghe,
Convex Optimization
.
Cambridge
University Press, Mar. 2004.
[13] S. Diamond and S. Boyd, “CVXPY: A Python-embedded modeling
language for convex optimization,”
J. Mach. Learn. Res.
, 2016.
[14] MOSEK ApS,
MOSEK Optimizer API for Python 9.0.105
, 2020.
[15] S.-J. Chung, S. Bandyopadhyay, I. Chang, and F. Y. Hadaegh, “Phase
synchronization control of complex networks of Lagrangian systems on
adaptive digraphs,”
Automatica
, vol. 49, no. 5, pp. 1148–1161, 2013.
[16] S. Singh, A. Majumdar, J. Slotine, and M. Pavone, “Robust online
motion planning via contraction theory and convex optimization,” in
IEEE Int. Conf. Robot. Automat.
, May 2017, pp. 5883–5890.
[17] M. Johansson and A. Rantzer, “Computation of piecewise quadratic
Lyapunov functions for hybrid systems,”
IEEE Trans. Autom. Control
,
vol. 43, no. 4, pp. 555–559, Apr. 1998.
[18] T. A. Johansen, “Computation of Lyapunov functions for smooth non-
linear systems using convex optimization,”
Automatica
, vol. 36, no. 11,
pp. 1617 – 1626, 2000.
[19] P. Giesl,
Construction of Global Lyapunov Functions Using Radial Basis
Functions
, ser. Lecture Notes in Mathematics.
Springer-Verlag Berlin
Heidelberg, 2007, vol. 1904.
[20] A. Papachristodoulou and S. Prajna, “On the construction of Lyapunov
functions using the sum of squares decomposition,” in
IEEE Conf. Decis.
Control
, vol. 3, Dec. 2002, pp. 3482–3487.
[21] E. M. Aylward, P. A. Parrilo, and J.-J. E. Slotine, “Stability and
robustness analysis of nonlinear systems via contraction metrics and
SOS programming,”
Automatica
, vol. 44, no. 8, pp. 2163 – 2170, 2008.
[22] I. R. Manchester and J. E. Slotine, “Control contraction metrics: Convex
and intrinsic criteria for nonlinear feedback design,”
IEEE Trans. Autom.
Control
, vol. 62, no. 6, pp. 3046–3053, Jun. 2017.
[23] W. Tan, “Nonlinear control analysis and synthesis using sum-of-squares
programming,” Ph.D. dissertation, Univ. California, Berkeley, 2006.
[24] J. R. Cloutier, “State-dependent Riccati equation techniques: An
overview,” in
Proc. Amer. Control Conf.
, vol. 2, 1997, pp. 932–936.
[25] T. C ̧ imen, “State-dependent Riccati equation (SDRE) control: A survey,”
IFAC Proc. Volumes
, vol. 41, no. 2, pp. 3761 – 3775, 2008.
[26] H. K. Khalil,
Nonlinear Systems
, 3rd ed. Prentice-Hall, 2002.
[27] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan,
Linear Matrix
Inequalities in System and Control Theory
, ser. Studies in Applied
Mathematics. Philadelphia, PA: SIAM, Jun. 1994, vol. 15.
[28] R. A. Horn and C. R. Johnson,
Matrix Analysis
, 2nd ed.
Cambridge
University Press, 2012.
[29] S. Bonnabel and J. Slotine, “A contraction theory-based analysis of
the stability of the deterministic extended Kalman filter,”
IEEE Trans.
Autom. Control
, vol. 60, no. 2, pp. 565–569, Feb. 2015.
[30] Y. Nakka
et al.
, “A six degree-of-freedom spacecraft dynamics simulator
for formation control research,” in
Astrodyn. Special. Conf.
, 2018.