Contracting Nonlinear Observers: Convex Optimization and
Learning from Data

arXiv:1711.08135v1 [cs.SY] 22 Nov 2017

Ian R. Manchester1
Abstract— A new approach to design of nonlinear observers
(state estimators) is proposed. The main idea is to (i) construct
a convex set of dynamical systems which are contracting
observers for a particular system, and (ii) optimize over this
set for one which minimizes a bound on state-estimation error
on a simulated noisy data set. We construct convex sets of
continuous-time and discrete-time observers, as well as contracting sampled-data observers for continuous-time systems.
Convex bounds for learning are constructed using Lagrangian
relaxation. The utility of the proposed methods are verified
using numerical simulation.

I. I NTRODUCTION
The problem of estimating the “internal state” (a.k.a.
hidden or latent variables) of a dynamical system is one of
the canonical problems in control engineering, and similar
problems are encountered in time-series prediction, machine
learning, signal processing, and may other fields. For linear
systems a comprehensive methodology is available that combines computational simplicity, stability, and certain types of
statistical or deterministic optimality (e.g. [1], [2]).
For non-linear systems, on the other hand, compromises
are necessary. In all but a few isolated cases there is no
finite-dimensional representation of the statistically-optimal
estimator, so one must sacrifice one or more of optimality,
global stability, or computational simplicity.
Extended Kalman filters retain the simple structure of the
linear case, but are suboptimal and stability is generally
difficult to establish [1]. More accurate but computationally
intensive approaches include particle filters a.k.a. sequential Monte Carlo methods [3], Moving-horizon estimation
[4], and set-valued state estimators [5]. Notably, particle
filters and set-valued estimators maintain high-dimensional
representations of a distribution (or set) of possible states,
rather than a simple point estimate, while moving-horizon
estimation maintains a long history of recorded inputs and
outputs. Therefore even for low order systems, the internal
state of the estimator is of high dimension.
The term “nonlinear observer” usually refers to a state
estimator that abandons any claim to statistical optimality,
but is a dynamical system of low dimension (equal or
similar dimension to the true system), which is simple to
implement, and for which global or semi-global stability can
be established analytically or computationally [6].
Over the years, many methods of observer design have
been developed, including those based on geometric methods
*This work was supported by the Australian Research Council
1 I. Manchester is with the Australian Centre for Field Robotics, School
of Aerospace, Mechanical and Mechatronic Engineering,, University of
Sydney, Australia. ian.manchester@sydney.edu.au

[7], high-gain designs [8], immersion and invariance [9], and
finding a nonlinear transformation (possibly with excessive
coordinates) to a stable linear observer [10].
Computational approaches based on convex optimization
(in particular, linear matrix inequalities and semidefinite
programming) have been used to construct observers for systems with nonlinearities characterized by monotonicity [11],
[12] and L2 gain (e.g. [13] and references therein). Most
approaches are based on the search for a quadratic Lyapunov
function, while in [14] sum-of-squares programming was
used to find observers for polynomial systems via certain
non-quadratic Lyapunov functions. Sum-of-squares relaxes
the search for sign-definite polynomials to a semidefinite
program [15], and will also be applied in the present paper.
Several papers have addressed observer design using contraction analysis, e.g. [16], [17], [18]. Indeed, observer design
was one of the original motivations for contraction analysis
as introduced in [16]. Contraction analysis is based on the
study of a nonlinear system by way of its differential dynamics (a.k.a. variational system) along solutions. Roughly
speaking, since the differential dynamics are linear timevarying (LTV), many techniques from linear systems theory
can be directly applied. A central result is that if all solutions
of a smooth nonlinear system are locally exponentially
stable in a common metric, then all solutions are globally
exponentially stable. Historically, basic convergence results
on contracting systems can be traced back to the results of
[19] in terms of Finsler metrics, further explored in [20],
while convex conditions for existence and robustness of limit
cycles were given in [21], and extensions to constructive
feedback design were given in [22].
In this paper we will construct convex sets of contracting observers for continuous and discrete-time nonlinear
systems, as well as sampled-data observers, i.e. discretetime observers for continuous-time systems. Sampled-data
observer design is complicated by the inability to exactly
solve a nonlinear differential equation over a sampling interval, and simple numerical integration techniques can cause
instability [23]. We will construct sampled-data observers via
trapezoidal approximation, and prove that they remain stable.
Having constructed a convex set of observers, we propose
optimizing over this set to find an observer that minimizes
a measure of state estimation error on a finite data set. That
is, our approach is based on learning from data, and is
similar to the prediction error methods in black-box nonlinear
system identification [24], [25] and recurrent neural networks
[26], and snapshot-based model reduction [?]. The central
difference is that in prediction-error method the objective

is to find a dynamical model of the system, whereas in
this paper we assume a dynamical model is known, and
the objective is to optimize a state estimator. This could
be considered a form of regularization of prediction error
methods: the identified model must behave as an observer
for a certain known system.
This paper builds upon a framework introduced for system
identification and model reduction in [27], [28]. In particular,
we make use of similar convex parameterizations of stable
systems and bounds on nonlinear system behaviour as those
papers. The novel contributions of this paper are: (i) we show
how these techniques can be applied to the observer design
problem, (ii) we provide new convex sets of models for CT
systems and exponentially contracting systems, (iii) we show
that trapezoidal integration of CT observers preserves stability, and (iv) we propose a method to optimize mean-square
error of observers based on simulated data. The benefits of
the proposed approach are illustrated using simulations.
II. P RELIMINARIES
A. Notation
In this paper, we will consider both continuous-time (CT)
and discrete-time (DT) systems and signals. Let Z+ and
R+ denote the non-negative integers and reals, respectively.
Let `n2 denote the set of square-summable signals Z+ →
Rn . Similarly, Ln2 denotes the subspace of square-integrable
signals R+ → Rn . We use subscripts for time indexing, i.e.
for a signal x, xt ∈ Rn denotes its value at time t ∈ Z+
or R+ . For a symmetric matrix, A > 0 (A ≥ 0) denotes
positive-definiteness (semidefiniteness), and A > B (A ≥ B)
means A − B > 0 (A − B ≥ 0). For a vector x, |x| is
the Euclidean norm, and |x|2P is shorthand for x0 P x where
P > 0. We will also abuse notation slightly by using the
shorthand |H|2P = H 0 P H even if H is matrix with more
than one column.
B. Observer Design
In this subsection we describe the deterministic nonlinear
observer design problem. In Section IV we will extend
consideration to stochastic models with disturbances and
measurement errors.
We consider a “true system” model with finite-dimensional
state x ∈ Rn , driven by a known input u ∈ Rm and
producing a measured output y ∈ Rp . The dynamic and
sensor models are given by the following equations
σxt = a(xt , ut ),

y = g(xt , ut ),

(1)

d
for CT models, and the
where the operator σ represents dt
shift operator xt 7→ xt+1 for DT models. In the CT case,
for simplicity of development we will assume existence and
uniqueness of solutions for all t ∈ [0, ∞).
In this paper, an observer for (1) is another dynamical
system with a state estimate as its internal state x̂ ∈ Rn ,
which takes the input and output of (1) as its own inputs:

σx̂t = φ(x̂t , ut , yt ),

(2)

While more general notions of an observer are possible with
different dimension to the true system [6], for this paper we
restrict attention to the simplest case when the observer’s
internal state is x̂. We will assume that the functions a, g, φ
are continuously differentiable in their arguments.
The observer design problem is to find a function φ in (2)
such that x̂t is uniquely-defined on t ∈ [0, ∞) (which we
call well-posedness), and guaranteed to converge (in some
sense) to xt as t → ∞. In this paper, we will focus on the
following forms of convergence:
Definition 1: System (2) is said to be an L2 observer for
(1) if for all x0 ∈ Rn , x̂0 ∈ Rn xt − x̂t ∈ Ln2 for CT or
`n2 for DT observers, and if x̂0 = x0 then x̂t = xt for all
t ≥ 0. System (2) is said to be an exponential observer with
rate λ > 0 if the stronger condition holds that there exists a
function b : Rn × Rn → R+ with b(x, x) = 0 for all x, such
that for all x0 , x̂0 ∈ Rn :
|xt − x̂t | ≤ e−λt b(x0 , x̂0 ).
(3)
C. Contraction Analysis
Contraction analysis [16] studies the convergence of solutions of a dynamical system to each other, rather than to a
particular pre-defined “nominal” solution as in standard Lyapunov analysis. Thus stability of a system can be established
independent of knowledge of a particular solution, a property
which is particularly well suited to observer design.
Contraction analysis is based on the study of the differential dynamics (a.k.a. variational, linearized, prolonged),
defined along solutions of a system of the form (2):
∂φ
δt .
σδt =
∂ x̂
A central result of [16] for the CT case is that if there exists
a uniformly bounded metric M (x) such that
∂φ 0
∂φ
M +M
≤ −2λM,
∂x
∂x
where Ṁ = ∂M
∂x φ, then the system is contracting with rate
λ. In the DT case we have the similar condition
∂φ 0
∂φ
M (xt+1 )
− M (xt ) ≤ (1 − e−λ )M (xt ).
∂x
∂x
D. Bounds and Identities for Quadratic Functions
Throughout the paper we will frequently use the following
simple property. Concave functions g(c, P ) = −c0 P −1 c with
c ∈ Rn , P = P 0 > 0 obey the upper bound:
Ṁ +

−c0 P −1 c ≤ b0 P b − 2b0 c

(4)

where the right-hand side is a convex function of c, P for
any fixed b ∈ Rn . Inequality (4) follows directly from the
expansion b0 P b − 2b0 c + c0 P −1 c = |c − P b|2P −1 ≥ 0. From
this expansion it is also clear that the bound (4) is tight if
c = P b.
We will also frequently use the following version of the
polarization identity, valid for Q = Q0 > 0, arbitrary
matrices E, F of appropriate dimension, and arbitrary scalar
h > 0, which follows from expansion of the right hand side:


E 0 QF + F 0 QE = h1 |E + h2 F |2Q − |E − h2 F |2Q . (5)

III. C ONVEX S ETS OF N ONLINEAR O BSERVERS
In this paper, a system set is a pair (φ, Θ) where φ :
Rn × Rm × Rp × Rq → Rn is a continuously differentiable
function, and Θ ⊂ Rq is a set of q-dimensional parameter
vectors. Associated with (φ, Θ) are state-space dynamical
systems of the form

In principle, our results apply to search over infinite dimensional space of functions e, f that are continuously
differentiable, but in practical implementations e and f will
be linearly parameterized by a finite-dimensional vector θ:
eθ (x̂) =

q
X

θi ei (x̂), fθ (x̂, u, y) =

i=0

σx̂t

= φ(x̂t , ut , yt , θ),

(6)

where x̂ ∈ Rn , u ∈ Rm , y ∈ Rp , θ ∈ Θ.
We define a contracting observer set for a system (1) as
any system set such that, for all θ ∈ Θ, the following two
conditions are satisfied:
1) Contraction: any pair trajectories x̂a , x̂b of (2) with the
same inputs (u, y) but with different initial conditions
x̂a0 , x̂b0 converge, in one of the following two senses
a) L2 contraction, i.e. x̂a − x̂b ∈ L2 (or `2 for DT
observers).
b) exponential contraction with rate λ > 0, i.e.
|x̂a0 − x̂b0 | ≤ e−λt b(x̂a0 , x̂b0 ).

(7)

function b : Rn × Rn → R+ with b(x, x) = 0.
2) Correctness: when initialised with x̂0 = x0 , the observer matches the true system, i.e. x̂t = xt for all
t ≥ 0. This is the case if and only if
a(x, u) = φ(x, u, g(x, u), θ) ∀x ∈ Rn , u ∈ Rm . (8)
In what follows, we will usually drop the θ argument, and
speak directly of searching over functions φ(x̂, u, y).
It is obvious that these two conditions are sufficient for
(2) to be an observer for (1) for any θ ∈ Θ: correctness
implies that the true state x is a particular solution of the
observer, while contraction implies that all solutions of the
observer converge to each other, hence all solutions of the
observer converge to the true state. This characterization
of a contracting observer is a special case of the idea of
virtual system used to study observers and more general
synchronization behaviors in [29].
A. Convex Sets of Contracting Continuous-Time Observers
In this section we will construct convex sets of contracting
nonlinear observers for CT systems of the form
ẋt = a(xt , ut ),

yt = g(xt , ut ).

(9)

We will construct our observers in the following implicit
representation:
d
e(x̂t ) = E(x̂t )x̂˙ t = f (x̂t , ut , yt ).
dt

(10)

where E(x̂) = ∂∂x̂ e(x̂). As noted in [27], [28] an implicit
system representation significantly expands the flexibility of
the system set while retaining convexity.
We will ensure that E(x̂) is invertible for all x̂, and hence
the observer can be rewritten in the explicit form
x̂˙ t = φ(x̂t , ut , yt ) = E(x̂t )−1 f (x̂t , ut , yt ).

(11)

q
X

θi fi (x̂, u, y), (12)

i=0

where θ ∈ Rq and each basis function ei : Rn →
Rn , fi : Rn ×Rm ×Rp → Rn , is continuously differentiable
in its arguments. In particular, if the basis functions are
polynomials, then sum-of-squares programming can be used
to guarantee global sign-definiteness [15].
To show contraction, we make use of the differential
dynamics of (10), defined along a particular solution x̂, u, y:
d
(E(x̂t )δt ) = F (x̂t , ut , yt )δt .
dt

(13)

where E(x̂) = ∂∂x̂ e(x̂) and F (x̂, ut , yt ) = ∂∂x̂ f (x̂, ut , yt ).
Following [27], we will consider contraction metrics of the
form V (x̂, δ) = δ 0 E(x̂)0 P −1 E(x̂)δ where P = P 0 > 0
is an auxiliary matrix variable. Contraction with respect
to this metric can be established by differential dissipation
inequalities:
d
V (x̂, δ) = δ 0 (E 0 P −1 F + F 0 P −1 E)δ ≤ −δ 0 Hδ (14)
dt
for all δ, x̂, u, y. In particular, L2 contraction follows from
any H > 0, while exponential contraction with rate λ follows
from the choice H = 2λE 0 P −1 E. The optimal bound
b(x̂a0 , x̂b0 ) in (7) is the Riemannian energy with respect to
V [22].
The correctness condition for observer (10) and system (9)
is the following linear (and hence convex) constraint:
E(x)a(x, u) = f (x, g(x, u), u) ∀x ∈ Rn , u ∈ Rm

(15)

With E invertible, this clearly implies (8) via (11).
So the objective is to find an observer (10) with E(x̂)
invertible for all x̂ that satisfies the contraction condition
(14) and the correctness condition (15). The difficulty is that
both invertibility of E(x̂) and (14) are nonconvex constraints
on θ. We first recall the following useful result [28, Thm. 1],
giving a convex constraint for invertibility:
∂e
Lemma 1: Let E(x) = ∂x
for continuously differentiable
function x 7→ e(x) mapping Rn → Rn . Suppose
E(x) + E(x)0 ≥ µI,

∀x ∈ Rn

(16)

for some µ > 0, then e is a bijection, and E is non-singular
for all x.
We now present two alternative choices of convex constraints guaranteeing the contraction condition (14).
1) CT1: We restrict the space of functions e to the class
of e(x) = P x. Then (14) reduces to the condition
F + F0 + H ≤ 0

(17)

with H > 0 for L2 contraction or H = 2λP for
exponential contraction with rate λ. This is jointly
convex in e, f, P and quasi-convex in λ.

2) CT2: Condition (14) is replaced with the stronger
condition
(E − h2 F ) + (E − h2 F )0 − P
−(E + h2 F )0 P −1 (E + h2 F ) − H ≥ 0.

(18)

for some scalar h > 0 and arbitrary matrix H > 0
for L2 contraction or H = 2λE 0 P E for exponential
contraction with rate λ.
Both of these conditions are jointly-convex in e, f, P , and
the exponential contraction conditions are quasi-convex in
λ > 0. Therefore the rate λ can be maximized via bisection
search with a convex feasibility problem at each step.
The following theorem summarizes the main results of this
section:
Theorem 1: Given a true system (9), a convex set of
contracting observers is given by the system set (10), (12),
with Θ defined by the convex constraints (16), (15), and
either (17) or (18), with H as specified for L2 or exponential
contraction.
Proof: As remarked above, well-posedness follows
from (16) via Lemma 1. The correctness condition (15) is
clearly equivalent to (8) via the explicit representation (11).
To prove that (17) implies contraction for CT1, note that
with e(x̂) = P x̂ the observer is
x̂˙ = P −1 f (x̂, u, y)
and contraction (14) follows directly from (17) since E = P .
To show that (18) implies contraction for CT2, we first
apply (5) to (14) gives the equivalent condition


δ 0 |E + h2 F |2P −1 − |E − h2 F |2P −1 δ ≤ −2hλδE 0 P −1 Eδ
Then applying (4) with c = (E − h2 F )δ and b = δ shows
that (18) guarantees (14).
Remark 1: By using different quantities for b in (4) one
obtains a family of convex sets guaranteeing contraction. For
example, an iterative refinement procedure could start with
CT1 or CT2 and take advantage of the fact that (4) is tight
when c = P b and set bk+1 = Pk−1 ck = Pk−1 (Ek − Fk ),
where k indexes iteration number, such that the convex bound
is tight at the result of the previous iteration.
B. Convex Sets of Contracting Discrete-Time Observers
This section contains parallel results for DT systems:
xt+1 = a(xt , ut , ),

yt = g(xt , ut ).

(19)

The observers we will consider are of the implicit form
eθ (x̂t+1 )

= fθ (x̂t , ut , yt ),

(20)

with e, f parameterized as in (12), but we will ensure that e
is a bijection, and hence the observer can be written in the
explicit form
x̂t+1 = φ(x̂t , ut , yt ) = e−1
θ (fθ (x̂t , ut , yt ))

(21)

Well-posedness, i.e. the fact that e is a bijection, will be
guaranteed by ensuring E+E 0 ≥ µI, and applying Lemma 1.

A solution can be generally be obtained rapidly by Newton’s
method.
The correctness condition in DT is
e(a(x, u)) = f (x, u, g(x, u)) ∀x ∈ Rn , u ∈ Rm .

(22)

This constraint is linear (and hence convex) in e and f for
any known functions a, g.
To study contraction, we make use of the differential
dynamics of (20), given by
Et+1 δt+1 = Ft δt

(23)

where we have used the shorthand Et := E(x̂t ) and Ft :=
F (x̂t , ut , yt ).
Similarly to the CT case, contraction for the DT system
will be evaluated via a metric of the form V (x̂, δ) =
δ 0 E(x̂)0 P −1 E(x̂)δ. The contraction condition is then
Vt+1 − Vt = δt0 (Ft0 P −1 Ft − Et0 P −1 Et )δt ≤ −δt0 Hδt (24)
again using shorthand Vt := V (x̂t , δt ), and H > 0 for `2
contraction, and H = (1 − e−λ )Et0 P −1 Et for exponential
convergence with rate λ. Condition (24) is not jointly convex
in e, f, P due to the concave term −Et0 P −1 Et .
Using (4) to convexify (24) with b = δt , C = Et δt we
obtain the following:
E(x̂) + E(x̂)0 − P − F (x̂, u, y)0 P −1 F (x̂, u, y) − H ≥ 0
(25)
for all x̂ ∈ Rn , u ∈ Rm , y ∈ Rp . To summarize the
construction in this section:
Theorem 2: Given a true system (19), a convex set of
contracting observers is given by the parameterization (20),
(12) where Θ is defined by the following constraints
1) Contraction: for some P > 0, (25) holds with H >
0 arbitrary for `2 contraction, and H = 2λE 0 P E for
exponential contraction with rate λ.
2) Correctness: (22) holds for all x, u.
which are convex in the unknowns e, f, P .
Note that well-posedness follows automatically from (25)
and Lemma 1.
C. Contracting Sampled-Data Observers for CT Systems
It is common in applications that a dynamical model is
known in continuous time (e.g. from physical laws), but an
observer is to be implemented in discrete time on a computer.
For nonlinear systems, the challenge is that exact solution
of CT differential equation is not available, so precise
sampling is generally not possible. Simple methods, such
as Euler integration, can destroy stability properties of CT
observers [23]. In this section we show that trapezoidal
integration:
h
zt+1 − zt ≈ (żt+1 − żt )
2
where h > 0 is the time-step, is particularly well-suited
to sampled-data implementation of contracting observers.
It has long been known that trapezoidal integration enjoys

1.5

favourable stability properties. In fact it is the most accurate linear multi-step method that guarantees stability when
integrating arbitrary linear systems [30].
Given a CT observer satisfying the conditions of Section
III-A, a sampled-data observer can be constructed as
h
h
f (x̂t+h , ut+h , yt+h ) = e(x̂t ) + f (x̂t , ut , yt ).
2
2
(26)
Note that the update from t → t + h depends on measured
data of u and y at both times t and t + h.
Theorem 3: For a CT system (9) and any CT observer
satisfying the conditions of Theorem 1, the sampled-data
observer given in (26) is well-posed, `2 stable, and “correct”
in the sense that it matches the trapezoidal integration of the
CT system
Proof: We sketch the proof due to space restrictions.
Well-posedness of the observer (26) is guaranteed by Lemma
1 if E(x̂)− h2 F (x̂, u, y)+E(x̂)0 − h2 F (x̂, u, y)0 ≥ µI. If CT1
is used, then by construction E + E 0 > 0 and F + F 0 ≤ 0
so clearly the observer is well-posed for any h. For CT2, by
inspection the only term on the left-hand-side of (18) that is
not negative-semidefinite is (E− h2 F )+(E− h2 F )0 . Therefore
since P > 0 this term is positive-definite, which is precisely
the condition required for well-posedness.
Both CT1 and CT2 imply (14). For the sampled data
observer consider the differential dynamics:

1

0.5

0

-0.5

e(x̂t+h ) −

(Et+1 − h2 Ft+h )δt+h = (Et + h2 Ft )δt
and metric
V (x̂, δ, t) = δ 0 (Et − h2 Ft )0 P −1 (Et − h2 Ft )δ
Then V (x̂t+h , δt+h , t + h) − V (x̂, δ, t) = δ(|Et − h2 Ft |2P −1 −
|Et + h2 Ft |2P −1 ). By the polarization identity and (14) this
difference is uniformly negative-definite, and hence the observer is `2 contracting
Note that exponential convergence with the same rate λ is
not guaranteed, since the metric used to evaluate contraction
is different.
D. Flexibility of the Observer Sets
To be able to find observers for a wide class of systems, it
is of course beneficial to have as flexible a set of observers
as possible. The main result regarding flexibility is the
following:
Theorem 4: Considering observers in the span of (12),
1) CT1, CT2, and DT contain all observers that can be
written in the form
σx̂t = Ao x̂t + q(ut , yt )
where Ao is a Hurwitz stable matrix for CT or Schur
stable matrix for DT.
2) CT1 and DT additionally contain all nonlinear observers
that are contracting with respect to a constant metric
V (x̂, δ) = δ 0 M δ where M = M 0 > 0.
we omit the proof due to space restrictions. Similar results
are proved in [27], [28].

-1

-1.5
-4

-3

-2

-1

0

1

2

3

4

Fig. 1. Graph of the function φ(x̂) = −(x̂ + x̂3 )/(0.1 + 3x̂2 ). The
positive slope around x̂ = ± 12 implies that the system x̂˙ = φ(x̂) is not
contracting with respect to a constant metric. However, it is in the set CT2.

This implies that observers exist in these sets for any
system of the form
σxt = Axt + q(ut , yt ),

yt = Cxt + r(ut )

where the pair (A, C) is detectable, including all full-order
observers for linear detectable systems.
From Theorem 1, it may appear that CT1 is more flexible
than CT2. In fact CT1 contains systems that CT2 does not,
and CT2 contains systems that CT1 does not. In particular,
considering polynomial bases for e, f , CT1 contains the
system:
x̂˙ = −x̂ − x̂3
with e(x̂) = x̂, f (x̂) = −x̂ − x̂3 , since F (x̂) = −1 − 3x̂2 is
uniformly negative-definite.
CT2 does not contain this system: It is clear from (18)
that (E − F ) + (E − F )0 must be larger in magnitude than
(E + F )0 P −1 (E + F ), which is impossible if E is constant
and F is unbounded in magnitude as x̂ 7→ ∞, as is the case
if f is a polynomial of degree > 1.
On the other hand, CT2 contains systems that CT1 does
not. In particular, the system
x̂ + x̂3
x̂˙ = φ(x̂) = −
0.1 + 3x̂2
from e(x̂) = 0.1x̂+ x̂3 , f (x̂) = −x̂− x̂3 . It can be seen from
the graph of φ in Fig. 1 that the slope is positive, and hence
∂φ 0
there is no constant P > 0 such that P ∂φ
∂ x̂ + ∂ x̂ P < 0.
However, by direct evaluation this system satisfies (14) with
P = 1, H = 0.1. However, it is contracting with respect to
the metric V (x̂, δ) = δ 0 E(x̂)0 P −1 E(x̂)δ = (0.1 + 3x̂2 )2 δ 2 .
IV. L EARNING FROM DATA
Within a set of stable observers, it is natural to try to
find the one that is “best” in some sense. A standard way
to formulate the problem is to specify a particular stochastic
model of the true system, e.g. in the DT case
xt+1 = a(xt , ut , wt ),

y = g(xt , ut , wt )

(27)

where ut is a known input, and wt represents measurement
and/or process noise, and is unknown but drawn from a
known probability distribution (we assume w = 0 corresponds to a “nominal” model). Initial conditions may be

known, or also drawn from a probability distribution. Then
a measure of state estimation error is chosen, e.g. the meansquare error:
T
1X
E[(xt − x̂t )0 (xt − x̂t )],
T →∞ T
t=0

J = lim

where E denotes the expectation operator (see, e.g., [1]).
Unfortunately there is no computationally tractable solution
to this problem for most nonlinear systems, so some form
of approximation is required.
We propose an approach based on “learning from data”.
The procedure we propose is outlined as follows:
1) Construct a “flexible” set of contracting nonlinear observers for the nominal (w = 0) model.
2) Draw one or more realisations of w and simulate
the stochastic model (27), collecting data sets z̃ =
{x̃t , ũt , ỹt |t = 0, 1, ..., T } as the realisations of x, u, y.
3) Optimize over the set of observers for one which
minimizes the empirical mean-square error:
J=

T
1X
|x̂t − x̃t |2
T t=0

(28)

˜ 0 = 0, where the equation errors t are defined as
with ∆
t := e(x̃t+1 , ut+1 , yt+1 ) − f (x̃t , ũt , ỹt )
Then the local mean-square error is
J 0 (θ, z̃) =

T
1X
|∆t |2 ,
T t=0

~ = ~(θ, z̃)
H(θ, z̃)∆
~ := [∆
˜0 ,∆
˜ 0 , ..., ∆
˜ 0 ]0 and
with the stacked vectors ∆
1
2
T
~(θ, z̃) := [00 , 01 , ...0T ]0 and H(θ, z̃) is a block lowerbidiagonal matrix that is affine in θ and is straightforward to
construct from (28).
Then a convex upper bound for J 0 is given by
~ 2 − 2∆
~ 0 (H(θ, z̃)∆
~ − ~(θ, z̃))}.
Jˆ0 (θ, z̃) := sup {|∆|
~ nT
∆∈l
T

It can be shown that for all systems satisfying the contraction
~ is finite, and can be written
constraint, the supremum over ∆
in explicit form:
Jˆ0 (θ, z̃) = ~(θ, z̃)0 H(θ, z)−1~(θ, z̃)

where x̂ are solutions of the observer simulated with
data ũt , ỹt , and an estimate of x̂0 . When multiple
realisations of w are sampled, J is the sum of the
corresponding mean-square errors.
If a method is available for simulating a stochastic CT model,
then the same approach can be applied to optimize sampleddata observers by subsampling at fixed rate.
Remark 2: In the linear case, the constraint that the
learned observer is an observer for the nominal system means
the resulting estimator is unbiased when w is zero mean. In
the nonlinear case this is not generally true, but it can be
considered as a form of regularization that introduces a bias
to ensure reliable behaviour on data not in the training set.
While the above recipe can be followed using generic
nonlinear programming methods to minimize J, this may
prove challenging since J is a highly non-convex function of
the observer parameters even in the linear case. We construct
a convex approximation based on Lagrangian relaxation of
linearized estimation error, similar to the method proposed
in [28] for system identification.
First we construct the linearization along the sampled “true
˜ t ≈ x̂t − x̃t
solution” x̃, by introducing the local deviation ∆
obeying the dynamics
˜ t + t .
E(x̃t+1 )∆t+1 = F (x̃t , ũt )∆

Roughly speaking, when the state estimation error is small,
˜ t ≈ x̂t − x̃t then J ≈ J 0 . A more formal derivation
i.e. ∆
for system identification problems can be found in [28, Sec.
V.B]. Note that if the observer dynamics are affine in x̂, then
J and J 0 are identical.
The problem of minimizing J 0 remains nonconvex, but a
semidefinite programming bound can be constructed using
Lagrangian relaxation [28]. This is most clearly expressed
be rewriting (28) in a “lifted” representation

(29)

(30)

where H = H + H 0 − I or, via Schur complement, in the
equivalent semidefinite programming representation:


s
~(θ, z̃)0
Jˆ0 (θ, z̃) = min s s.t.
≥ 0.
(31)
~(θ, z̃) H(θ, z)
Algorithms have been developed for efficient solution of
problems with this structure [31], [32]. Both have computational complexity that is linear in the data length T , rather
than cubic for a generic semidefinite programming solver.
The main result regarding Jˆ0 is the following:
Theorem 5: For all P = P 0 > 0 and signal data z̃, for
every θ ∈ Θ, J 0 (θ, z̃) ≤ Jˆ0 (θ, z̃) < ∞. Furthermore, if the
data set z̃ is generated with w = 0, then the optimal values
satisfy J 0,? (z̃) = Jˆ0,? (z̃) = 0.
We omit the details of the proof because of space restrictions,
since it is similar to [28, Theorem 6].
V. E XAMPLES AND D ISCUSSION
In this section we illustrate the method with some numerical simulations. The first example system we consider was
previously studied in [12]:


x2
ẋ = a(x) =  x2 − 31 x32 − x2 x23  , y = g(x) = x1
x2 − x3 − 31 x33 − x3 x22
The solution given in [12] relied on a certain ingenuity
on the part of the designer, recognising the system can be
decomposed into a linear part and a monotonic nonlinearity.
In comparison, our method is quite “plug and play”: we apply
Theorem 1 using CT1 as a contracting set and f the set of
degree-three polynomials in x̂, y, the contraction constraint
is imposed via sum-of-squares [15], and an observer is found
in 1.41 seconds on 3.1 GHz Intel Core i7 with 16GB RAM.

4

4

3

3

2

2

1

States

States

1
0

-1

0

-1
-2

-2

-3

-3

-4

-5
0

0.5

1

1.5

2

2.5

3

-4
0

0.5

1

1.5

Time (s)

Fig. 2. Comparison of true states (solid lines), CT observer (dashed lines),
and sampled-data observer (circles) with sampling interval h = 0.1s.

2

2.5

3

Time (s)

Fig. 3. Comparison of true states (solid lines), CT observer (dashed lines),
and sampled-data observer (circles) with sampling interval h = 0.3s.
1.1

A bisection search over λ reveals that λ = 1 is the best
rate of exponential convergence in this observer set, and with
λ = 1 we computed the observer that minimized the `1 norm
of the coefficients of f , to encourage sparsity of coefficients.
All computations were done using Yalmip [33] and Mosek.
Fig. 2 shows the evolution of the three states (solid lines)
from initial conditions x0 = [−3, 4, 1]0 and the observer
states (dashed lines) initialized at x̂0 = [0, 0, 0]0 . Also shown
are the states of the sampled-data observer (circles) proposed
in Section III-C with a sampling interval h = 0.1s. It can be
seen that the sampled-data observer is a good approximation
to the CT observer, and converges rapidly to the true states.
In Fig. 3 we show corresponding results with a longer
sampling interval of h = 0.3s. It can now be seen that the
behaviour of the CT and sampled-data observers are quite
different in the initial transient, nevertheless both converge
to the true states.
We next tested the ability of the learning approach in
Section IV to find an observer with good performance on
noisy data. To generate training data we simulated the CT
system over 5 seconds, sampled the state x̃ and output with
sampling interval h = 0.1s, added Gaussian white noise of
variance 0.01 to the output to generate ỹ, and then minimized
Jˆ0 over the set CT1 with H = 10−8 I.
We tested the resulting observer against the previouslydesigned observer (optimal with respect to λ) on 50 realisations of validation data, generated in the same way. A
boxplot of the resulting mean-square errors is shown in Fig.
4. For this example, the learned observer has around half
the median error and one third the worst-case error of the
observer optimized for λ.
To benchmark the proposed learning method against a
“gold standard”, we now turn to observer design for linear
discrete-time systems and compare to the Kalman filter. We
compare performance on randomly sampled systems, using
different amounts of training data, as follows:
1) For each data length T ∈ {100, 250, 500, 1000, 2000}
do the following

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
Optimal Exponential Rate

Optimized via Lagrangian Relaxation

Fig. 4. State estimation error on noisy validation data, comparing the
observer which is optimal in terms of exponential rate to the one which is
optimized for noise performance via Lagrangian relaxation.

2) Sample 20 random 4th-order systems using Matlab
drss command, eliminating systems with pure integrators. For each system compute the steady-state Kalman
filter and do the following:
3) Generate a random training data set of length T and
train a DT observer using the method of Sec. IV.
4) Sample 25 random validation data sets of length and
compute the relative error for each, defined as a percentage:
Jrel = 100

PT

t=0

P
o
2
|x̂o
−x̃t |2 − T
tP
t=0 |x̂t −x̃t |
T
k −x̃ |2
|x̂
t
t
t=0

where x̂o is the state estimate from our learned observer,
and x̂k is the state estimate from the Kalman filter.
Boxplots of the computed Jrel for all realisations are
shown in Fig. 5, and the median values of Jrel for each
training-data length are collected in Table I. It can be
seen that with just a few hundred training samples, the
performance of the learned estimator comes within around
2% of the Kalman filter, and on some data sets outperforms
the Kalman filter, indicated by negative quantities in Fig 5.
VI. C ONCLUSIONS
In this paper, we have introduced methods for constructing
convex sets of contracting CT, DT, and sampled-data observers for a given nonlinear system. We propose optimizing

80

Relative variance (%)

60

40

20

0

-20

50

100

200

300

400

Number of training samples

Fig. 5. Comparison to the Kalman filter: boxplot of Jrel vs number of
samples used for training
TABLE I
M EDIANS OF Jrel FOR EACH TRAINING - DATA LENGTH
Samples (N)
Median Jrel (%)

50
8.4336

100
4.2223

200
3.4214

300
1.4955

400
1.3731

over these sets to based on sampled simulation data to “learn”
and observer with good performance on noisy data.
While we have focused on the case that the observer
state is x̂, the proposed methods can be adapted to reducedorder and excessive-order observers. Also important is the
case of “reduced-complexity” nonlinear observers, for which
the correctness condition is infeasible, i.e. when the true
system dynamics are not in the span of (12), one can relax
correctness to bounded error in (8) subsets of Rn × Rm .
Investigation of these approaches is underway.
R EFERENCES
[1] B. D. O. Anderson and J. B. Moore, Optimal Filtering. Prentice Hall,
1979.
[2] I. R. Petersen and A. V. Savkin, Robust Kalman Filtering for Signals
and Systems with Large Uncertainties. Springer Science & Business
Media, 1999.
[3] A. Doucet, N. De Freitas, and N. Gordon, Sequential Monte Carlo
Methods in Practice. Springer Verlag, 2001.
[4] C. V. Rao, J. B. Rawlings, and D. Q. Mayne, “Constrained state
estimation for nonlinear discrete-time systems: Stability and moving
horizon approximations,” IEEE transactions on automatic control,
vol. 48, no. 2, pp. 246–258, 2003.
[5] J. S. Shamma and K.-Y. Tu, “Approximate set-valued observers for
nonlinear systems,” IEEE Transactions on Automatic Control, vol. 42,
no. 5, pp. 648–658, May 1997.
[6] L. Praly, “Observers for Nonlinear Systems,” Encyclopedia of Systems
and Control, pp. 935–943, 2015.
[7] A. J. Krener and A. Isidori, “Linearization by output injection and
nonlinear observers,” Systems & Control Letters, vol. 3, no. 1, pp.
47–52, Jun. 1983.
[8] H. K. Khalil and L. Praly, “High-gain observers in nonlinear feedback
control,” International Journal of Robust and Nonlinear Control,
vol. 24, no. 6, pp. 993–1015, Apr. 2014.
[9] D. Karagiannis, D. Carnevale, and A. Astolfi, “Invariant Manifold
Based Reduced-Order Observer Design for Nonlinear Systems,” IEEE
Transactions on Automatic Control, vol. 53, no. 11, pp. 2602–2614,
Dec. 2008.
[10] V. Andrieu and L. Praly, “On the Existence of a Kazantzis–
Kravaris/Luenberger Observer,” SIAM Journal on Control and Optimization, vol. 45, no. 2, pp. 432–456, Jan. 2006.

[11] M. Arcak and P. Kokotović, “Nonlinear observers: A circle criterion
design and robustness analysis,” Automatica, vol. 37, no. 12, pp. 1923–
1930, Dec. 2001.
[12] X. Fan and M. Arcak, “Observer design for systems with multivariable
monotone nonlinearities,” Systems & Control Letters, vol. 50, no. 4,
pp. 319–330, Nov. 2003.
[13] D. Coutinho, C. de Souza, K. Barbosa, and A. Trofino, “Robust Linear
$H {\infty}$ Filter Design for a Class of Uncertain Nonlinear Systems: An LMI Approach,” SIAM Journal on Control and Optimization,
vol. 48, no. 3, pp. 1452–1472, Jan. 2009.
[14] C. Ebenbauer, J. Renz, and F. Allgower, “Polynomial Feedback
and Observer Design using Nonquadratic Lyapunov Functions,” in
Proceedings of the 44th IEEE Conference on Decision and Control,
Dec. 2005, pp. 7587–7592.
[15] P. A. Parrilo, “Semidefinite programming relaxations for semialgebraic
problems,” Mathematical Programming, vol. 96, no. 2, pp. 293–320,
2003.
[16] W. Lohmiller and J.-J. E. Slotine, “On Contraction Analysis for Nonlinear Systems,” Automatica, vol. 34, no. 6, pp. 683–696, Jun. 1998.
[17] R. Sanfelice and L. Praly, “Convergence of Nonlinear Observers on
With a Riemannian Metric (Part I),” IEEE Transactions on Automatic
Control, vol. 57, no. 7, pp. 1709–1722, Jul. 2012, 00014.
[18] A. P. Dani, S. J. Chung, and S. Hutchinson, “Observer Design
for Stochastic Nonlinear Systems via Contraction-Based Incremental
Stability,” IEEE Transactions on Automatic Control, vol. 60, no. 3,
pp. 700–714, Mar. 2015.
[19] D. Lewis, “Metric properties of differential equations,” American
Journal of Mathematics, pp. 294–312, 1949.
[20] F. Forni and R. Sepulchre, “A Differential Lyapunov Framework
for Contraction Analysis,” Automatic Control, IEEE Transactions on,
vol. 59, no. 3, pp. 614–628, 2014.
[21] I. R. Manchester and J.-J. E. Slotine, “Transverse Contraction criteria
for existence, stability, and robustness of a limit cycle,” Systems &
Control Letters, vol. 62, pp. 32–38, 2013.
[22] I. R. Manchester and J. J. E. Slotine, “Control Contraction Metrics:
Convex and Intrinsic Criteria for Nonlinear Feedback Design,” IEEE
Transactions on Automatic Control, vol. 62, no. 6, pp. 3046–3053,
Jun. 2017.
[23] M. Arcak and D. Nešić, “A framework for nonlinear sampled-data
observer design via approximate discrete-time models and emulation,”
Automatica, vol. 40, no. 11, pp. 1931–1938, Nov. 2004.
[24] J. Sjöberg, Q. Zhang, L. Ljung, A. Benveniste, B. Delyon, P.-Y.
Glorennec, H. Hjalmarsson, and A. Juditsky, “Nonlinear black-box
modeling in system identification: A unified overview,” Automatica,
vol. 31, no. 12, pp. 1691–1724, 1995.
[25] L. Ljung, System Identification: Theory for the User, 3rd ed. Englewood Cliffs, New Jersey, USA: Prentice Hall, 1999.
[26] R. Pascanu, T. Mikolov, and Y. Bengio, “On the difficulty of training
recurrent neural networks,” in International Conference on Machine
Learning, 2013, pp. 1310–1318.
[27] M. M. Tobenkin, I. R. Manchester, J. Wang, A. Megretski, and
R. Tedrake, “Convex optimization in identification of stable non-linear
state space models,” in 49th IEEE Conference on Decision and Control
(CDC). IEEE, 2010.
[28] M. M. Tobenkin, I. R. Manchester, and A. Megretski, “Convex
Parameterizations and Fidelity Bounds for Nonlinear Identification and
Reduced-Order Modelling,” IEEE Transactions on Automatic Control,
vol. 62, no. 7, pp. 3679–3686, Jul. 2017.
[29] W. Wang and J.-J. E. Slotine, “On partial contraction analysis for
coupled nonlinear oscillators,” Biological Cybernetics, vol. 92, no. 1,
pp. 38–53, Dec. 2004.
[30] G. G. Dahlquist, “A special stability problem for linear multistep
methods,” BIT Numerical Mathematics, vol. 3, no. 1, pp. 27–43, Mar.
1963.
[31] J. Umenberger and I. R. Manchester, “Specialized Algorithm for
Identification of Stable Linear Systems using Lagrangian Relaxation,”
in Proc. 2016 American Control Conference, Boston, MA, 2016.
[32] M. S. Andersen, J. Dahl, and L. Vandenberghe, “Implementation of nonsymmetric interior-point methods for linear optimization
over sparse matrix cones,” Mathematical Programming Computation,
vol. 2, no. 3-4, pp. 167–201, Dec. 2010.
[33] J. Lofberg, “Yalmip : A toolbox for modeling and optimization in
MATLAB,” in Proceedings of the CACSD Conference, Taipei, Taiwan,
2004.

