Convex Parameterizations and Fidelity Bounds for
Nonlinear Identification and Reduced-Order Modelling

arXiv:1701.06652v1 [] 23 Jan 2017

Mark M. Tobenkin

Ian R. Manchester

Abstract—Model instability and poor prediction of long-term
behavior are common problems when modeling dynamical systems using nonlinear “black-box” techniques. Direct optimization
of the long-term predictions, often called simulation error minimization, leads to optimization problems that are generally nonconvex in the model parameters and suffer from multiple local
minima. In this work we present methods which address these
problems through convex optimization, based on Lagrangian relaxation, dissipation inequalities, contraction theory, and semidefinite programming. We demonstrate the proposed methods with
a model order reduction task for electronic circuit design and
the identification of a pneumatic actuator from experiment.

I. I NTRODUCTION
Building approximate models of dynamical systems from
data is a ubiquitous task in the sciences and engineering.
Reduced order modelling and system identification are two
well-known cases, which differ mainly in the source of the
data. In reduced order modelling the data come from simulation “snapshots” of a complex model1 , which is usually
based on physical principles (e.g. [2], [3], [4]). The central
problem is finding a model that simulates accurately despite
unmodelled dynamics. In system identification the data come
from experimental recordings [5], [6], [7], [8], and this brings
the additional problems of measurement noise and exogenous
disturbances. When some known structure can be incorporated
into a model it is usually preferable to do so, but black-box
modeling plays an important role when first-principles models
are either weakly identifiable, too complicated for the eventual
application, or simply unavailable.
The basic objective we consider is to obtain a dynamical
system that accurately reproduces the observed relationship
between some given input and output data, and can reliably
generalize to other inputs. A complete solution using optimization methods requires a specification of a model class, a
choice of objective function, and ultimately an algorithm for
minimizing the objective function over the model class.
The model class we propose consists of discrete-time nonlinear state space models. Such models have previously been
studied in the context of model reduction (e.g. [3]) and system
identification (e.g. [9]). It includes all linear time-invariant
M. Tobenkin and A. Megretski are with the Department of Electrical
Engineering and Computer Science, Massachusetts Institute of Technology,
Cambridge, MA, 02139. Email: {mmt, ameg}@mit.edu
I. Manchester is with the Australian Centre for Field Robotics and School of
Aerospace, Mechanical and Mechatronic Engineering, University of Sydney,
NSW 2006, Australia. Email: ian.manchester@sydney.edu.au.
This was supported by National Science Foundation Grant No. 0835947,
the LANL Information Science and Technology Institute, and the Australian
Research Council (DP130100551).
1 We note that another class of model reduction algorithms work directly
from the defining equations of the high-order model rather than simulation
snapshots, e.g. balanced truncation and Hankel-optimal reduction [1].

Alexandre Megretski

(LTI) systems, systems defined by LTI internal dynamics
with nonlinear input/output mappings (e.g. Wiener and Hammerstein models [10]), and includes nonlinear autoregressive
models (e.g. [5], [11]) by choosing the state to be a truncated
history of outputs.
An alternative is to model the current output as a static function of a finite history of inputs, e.g. Nonlinear Finite-ImpulseResponse (NFIR) and Volterra models [5], [12]. However,
such models are inefficient at representing inherently dynamic
behaviour such as resonance, which is naturally expressed by
models with internal feedback.
The inclusion of feedback in a model means that stability
is no longer guaranteed, and in many applications an unstable
model would be unacceptable. The ability of a model to
“generalize” to inputs not in the training data depends critically
on model stability. Ensuring stability of models with feedback
is a significant problem and has been studied for many years
for linear system identification [8]. It is known that stability of
linear autoregressive models can be guaranteed if the input is
generated by an autoregressive process and the data sequence
is sufficiently long [13]. Stability of linear state space models
generated by subspace identification was addressed by padding
a shifted state matrix with zeros in [14], and by regularization
in [15]. A linear matrix inequality (LMI) parametrization of
stable linear models was given in [16].
The challenge of ensuring stability of nonlinear models with
feedback is discussed in [7], [8], [11]. To the authors’ knowledge, there are few prior published methods that guarantee
stable models regardless of the source of the data. Worthy
of note are the recent results of [17] giving small-gain and
passivity-like stability guarantees for reduction of a high-order
linear system in feedback with a convergent nonlinear system.
The objective function we focus on in this paper is simulation error, a.k.a. output error. This function is in general highly
non-convex for models with feedback, and it is challenging to
compute the global minimum [8]. Indeed, one can find poor
local minima for the class of second-order linear models [18].
A computationally simpler alternative is to minimize the
equation error, i.e. the one-step ahead prediction error, rather
than the long-term simulation error. Nonlinear ARX and
subspace identification, for example, use these criteria [5],
[6]. In many cases, the minimization amounts to basic least
squares and can be solved very efficiently. The disadvantage
is that there is in general no connection between the equation
error and the simulation error, since an unstable (or weakly
stable) model might be found for which small equation errors
accumulate under recursive simulation.
The main contributions of this paper are a class of nonlinear state-space models that is convex in its parameters and
guarantees stability, and a family of convex upper bounds

for minimization of simulation error. The combination of
these allows problems in reduced order modelling and identification to be approached using semidefinite programming.
Our approach builds upon the ideas proposed in [19] for
nonlinear autoregressive models, which was extended and
applied to electronics model reduction in [20]. The model set
and some fidelity bounds in this paper were reported without
proof in [21]. The present paper significantly extends these
results by introducing more accurate fidelity bounds based on
Lagrangian relaxation, complete proofs of technical results,
and new experimental results. We do not discuss the effects
of measurement noise, but we note that a modified version of
the proposed method gives consistent estimates when repeated
measurements are available [22], [23], [24].
II. P ROBLEM S ETUP & P RELIMINARIES
Spaces of vector signals, i.e. mappings {0, 1, ..., T } → Rn
are denoted lTn . For a symmetric matrix A = A0 , A ≥ 0
and A > 0 mean that A is positive semidefinite and positivedefinite, respectively. Inequalities A ≥ B, A > B mean
A − B ≥ 0, A − B > 0, respectively. For symmetric
positive-semidefinite Q the notation |v|2Q is shorthand for
v 0 Qv. Throughout the paper, signals with a tilde, e.g. ỹ(t),
refer to recorded data, whereas signals without, e.g. y(t) refer
to the corresponding quantity generated by an identified model.
Assume that training data are given in the form of finite
sequences of inputs and outputs: ũ = {ũ(0), ũ(1), ..., ũ(T )} ∈
lTm , ỹ = {ỹ(0), ỹ(1), ..., ỹ(T )} ∈ lTp . We use the notation z̃ :=
{ũ, ỹ} to represent the “available data”.
In this paper, a model class is a triple M = (a, g, Θ) where
a : Rn × Rm × Rq → Rn and g : Rn × Rm × Rq → Rp are
continuously differentiable functions, and Θ ⊂ Rq is a set of
q-dimensional parameter vectors. A model class (a, g, Θ) is
associated with state-space dynamical systems of the form

n

x(t + 1)

=

a(x(t), u(t), θ),

(1)

y(t)

=

g(x(t), u(t), θ),

(2)

m

p

where x ∈ R , u ∈ R , y ∈ R , θ ∈ Θ. It is desirable that
identified models generalize to inputs not in the training data,
it is appropriate to impose a strong notion of stability that
which ensures predictable response to a wide variety of inputs
independent of model initial conditions, such as the following:
Definition 1: A system (1), (2) is incrementally `2 stable if
for every pair of initial conditions x1 (0), x2 (0) and every input
sequence u(t), t = 1, 2, ..., the states x1 (t), x2 (t) and outputs
yP
1 (t) = g(x1 (t), u(t)) and yP
2 (t) = g(x2 (t), u(t)) satisfy
∞
∞
2
2
|x
(t)
−
x
(t)|
<
∞,
1
2
t=0 |y1 (t) − y2 (t)| < ∞.
t=0
In addition to a model set, one must define a measure
of model fidelity. There are many possibilities, but perhaps
the most natural and straightforward measure of quality is
simulation error, also known as output error:
Jse (θ, x, z̃) :=

T
X
t=0

|y(t) − ỹ(t)|2

(3)

which is to be minimized over parameters θ and states x for
given data z̃, subject to dynamic constraints (1), (2).

The major difficulties we seek to address in this paper are:
(a) Obtaining a useful set of stable models. For all non-trivial
parameterizations of functions a(x(t), u(t), θ), an explicit
and tractable characterization of the set of all θ defining
a stable model is unknown.
(b) Efficient global minimization of Jse . Even if the cost
function J(θ, x, z̃) is convex in θ, the dynamic constraints
(1) generally define a nonconvex set: the simulation
output is generated by the recursion
y(t) = g(a(a(...a(x(0), ũ(0), θ), ..., ũ(t − 1), θ), ũ(t), θ),
so even if a and g are linear in the parameter θ, Jse is a
highly nonlinear function of the model
In this paper, we construct a family of computationally
tractable approximations to the simulation error minimization problem. More specifically, we formulate model classes
M̂ = (a, g, Θ̂) and optimization problems of the form:
ˆ z̃), with the following properties:
minθ∈Θ̂ J(θ,
(a) M̂ is a finite-dimensional convex parametrization of
models satisfying the following inclusions: LT I ⊆ M̂ ⊆
M`2 where LT I refers to the set of all stable linear
time invariant models of appropriate dimension, and M`2
is the set of all (possibly nonlinear) models that are
incrementally `2 stable.
ˆ z̃) ≥
(b) Jˆ is a convex function of θ such that: (i) J(θ,
ˆ
Jse (θ, x, z̃) for any feasible x; (ii) J(θ, z̃) always has
ˆ z̃) = 0 when z̃ are
a finite minimum over θ; (iii) J(θ,
generated by simulating the model represented by θ. In
other words, the method perfectly recovers a “true” model
if one exists in the model class.
A key tool we use is the Lagrangian relaxation, or SProcedure, which is a method for approximating optimization
problems with “difficult” nonconvex constraints by “easier”
unconstrained problems with multipliers. Versions of this
method have been used extensively in robust control and
combinatorial optimization [25], [26]. In the context of this
paper the idea is to approximate problems of the form
J ? = min J(θ, x), subject to h(θ, x) = 0
x,θ

(4)

where h(θ, x) is a vector of constraint functions. In the
context of simulation-error minimization, J(θ, x) represents
the measure of model fit, e.g. (3), and h(θ, x) represents the
dynamical constraints, e.g. (1).
Suppose that J is a convex function of θ and h is an affine
function of θ, but both are nonconvex functions of θ and
x jointly. In our context, Lagrangian relaxation refers to the
approximation of the above problem by the following:
Jˆλ? = min Jˆλ (θ),
θ

Jˆλ (θ) := sup{J(θ, x) + λ(x)0 h(θ, x)}
x

(5)
where λ is some vector function analogous to a Lagrange
multiplier in nonlinear programming. The function Jˆλ (θ) is
clearly convex, since it is the supremum of a family of convex
functions [27]. It is also an upper bound for J(θ, x) for any
feasible x: if h(θ, x) = 0 in (5), then Jˆλ (θ) = J(θ, x), so the
supremum over x must be no smaller.

It is apparent that if the Lagrangian relaxation (5) is applied
using two different constraint functions that nevertheless have
the same feasible set: h1 (θ, x) = 0 ⇔ h2 (θ, x) = 0,
but which differ for other values of θ, x, then the resulting
values of Jˆλ (θ) can differ. Consequently, performance of
the Lagrangian relaxation may be improved by allowing a
redundant parametrization of dynamical systems.
We will frequently deal with convex functions of the form
f (c, P ) = c0 P −1 c, where c and P are decision variables,
c is a vector and P is a symmetric positive-definite matrix
of appropriate size. These functions can be represented by a
slack variable and an LMI using the Schur complement [27].
Concave functions g(c, P ) = −c0 P −1 c with P > 0 obey the
following upper bound:
− c0 P −1 c ≤ b0 P b − 2b0 c

(6)

where the right-hand side is a convex function of c, P for
any fixed b. Inequality (6) 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 (6) is tight if c = P b.
III. C ONVEX PARAMETERIZATION OF S TABLE M ODELS
The models we propose are constructed from three linearly
parametrized vector functions:
eθ (x) =

q
X

θi ei (x),

fθ (x, u) =

i=0

gθ (x, u) =

q
X

q
X

θi fi (x, u),

i=0

θi gi (x, u),

(7)

i=0

where θ ∈ Rq and each basis function ei : Rn → Rn , fi :
Rn × Rm → Rn , gi : Rn × Rm → Rp is continuously
differentiable in its arguments. These functions are associated
with a parametrization of implicit state space dynamic models
eθ (x(t + 1))
y(t)

=

fθ (x(t), u(t)),

(8)

= gθ (x(t), u(t)).

(9)

A particular choice of state dimension and bases for e, f, g
defines a parametrization of models of the form (1), (2) given
by a(x, u, θ) = e−1
θ (fθ (x, u)) and g(x, u, θ) = gθ (x, u), as
long as it can be guaranteed that eθ (x) is a bijection for all θ.
We denote by E, F, G the span of the functions (7) over
θ ∈ Rq . We assume that E, F, G contain all linear functions
x 7→ Ex, (x, u) 7→ Ax + Bu, (x, u) 7→ Cx + Du, respectively,
where A, B, C, D, E are matrices of appropriate dimensions,
and that F is closed under multiplication by a non-singular
n × n matrix, i.e. if f (·, ·) ∈ F then M f (·, ·) ∈ F for
all M ∈ Rn×n with det(M ) 6= 0. Simple choices for
E, F, G satisfying these assumptions are vectors of multivariate
polynomials or trigonometric polynomials up to some fixed
degree, parametrized by their coefficients. For the sake of
brevity, we will usually drop the subscript θ from the functions
e, f and g and speak of searching directly for e, f and g.
To ensure model stability, we propose the following constraint, requiring an auxiliary matrix variable P = P 0 > 0:
|f (x1 , u) − f (x2 , u)|2(P −1 ) − 2(x1 − x2 )0 (e(x1 ) − e(x2 ))
+ |x1 − x2 |2P +µI + |g(x1 , u) − g(x2 , u)|2 ≤ 0

(10)

for all x1 , x2 ∈ Rn , u ∈ Rm for some fixed µ > 0. Note that
this constraint is convex with respect to e, f, g, P and hence
convex with respect to a parameter vector θ when e, f, g, P
are affinely parametrized.
Definition 2: For a given choice of state dimension and
basis functions ei , fi , gi , let Θ be the set of θ ∈ Rq for
which the model (8), (9) satisfies the inequality (10) for
some P = P 0 > 0. Additionally, let M̂ be the model class
(e−1
θ (fθ (x, u)), gθ (x, u), Θ).
Model well-posedness, i.e. existence of a unique solution,
is guaranteed by the following theorem:
Theorem 1: Given a model (8), (9), if (10) holds for all
x1 , x2 , u, then e is a bijection Rn ↔ Rn .
Proof. We need to show that for every z ∈ Rn , the equation
e(x) = z has a unique solution. Firstly, note that (10) implies
2(x1 − x2 )0 (e(x1 ) − e(x2 )) ≥ µ|x1 − x2 |2 . It is known (see,
e.g., [28, Th. 18.15]) that if T : Rn → Rn is a nonlinear
operator satisfying (x1 − x2 )0 (T (x1 ) − T (x2 )) ≥ c|x1 − x2 |2
for all x1 , x2 ∈ Rn with c > 0, then the nonlinear equation
T (x) = 0 has a unique solution. Take T (x) = e(x) − z and
c = µ/2 and the theorem is proved.
For simulation of the model, a solution to the equation
e(x) = z can be found by application of the Newton method.
Convergence can be guaranteed by guarding with the ellipsoid
method [24, Sec 2.2.2].
Theorem 2: For any model of the form (8), (9), the condition (10) implies global incremental `2 stability.
Proof. Consider the incremental Lyapunov function candidate
V (x1 , x2 ) = |e(x1 ) − e(x2 )|2(P −1 ) with P = P 0 > 0. The
inequality
−|e(x1 ) − e(x2 )|2(P −1 ) + |f (x1 , u) − f (x2 , u)|2(P −1 )

+µ|x1 − x2 |2 + |g(x1 , u) − g(x2 , u)|2 ≤ 0
(11)

for all x1 , x2 ∈ Rn , u ∈ Rm is a sufficient condition for
global incremental `2 stability. Now apply the bound (6), with
c = e(x1 ) − e(x2 ) and b = x1 − x2 , giving

−|e(x1 )−e(x2 )|2(P −1 ) ≤ |x1 −x2 |2P −2(x1 −x2 )0 (e(x1 )−e(x2 )).

and substitute the right-hand-side of this inequality for
−|e(x1 )−e(x2 )|2(P −1 ) in (11). It is then clear that (10) implies
(11) and hence stability of the model.
Remark 1: It is straightforward to impose other behavioral
constraints, e.g. passivity, by the same convexification procedure and appropriate choice of supply rate.
Definition 3: A system x(t + 1) = a(x(t), u(t)), y(t) =
g(x(t), u(t)), defined by functions a : Rn × Rm 7→ Rn and
g : Rn × Rm 7→ Rp is called quadratically incrementally `2
stable, or more briefly: quadratically stable, if there exists an
M ∈ Rn×n with M = M 0 > 0 such that the condition

|a(x1 , u)−a(x2 , u)|2M −|x1 −x2 |2M ≤ −|g(x1 , u)−g(x2 , u)|2 ,
(12)
holds for all x1 , x2 ∈ Rn and u ∈ Rm . I.e. incremental
stability can be verified by a quadratic Lyapunov function.
Theorem 3: The model set M̂ includes all quadratically
incrementally `2 stable models of the form (1), (2) for which
a(·, ·) ∈ F and g(·, ·) ∈ G.

Proof. Suppose (12) holds for some M = M 0 > 0. Choose
e(x) = M x, f (x, u) = M a(x, u), and P = M , i.e. c = P b in
(6), and leave g(x, u) as is. Clearly this model is equivalent to
(1), (2), and by assumption that a(·, ·) ∈ F and that F is closed
under multiplication by M , f (·, ·) ∈ F. Direct substitution
into (10) reveals that if (12) holds then so does (10).
Remark 2: All stable linear systems are automatically stable in the above sense, as are models with stable linear
dynamics but nonlinear input and/or output maps, e.g. Wiener
and Hammerstein models, as long as g(x, u) is globally
Lipschitz in x.
IV. C ONVEX U PPER B OUNDS FOR S IMULATION E RROR
Let us suppose there is method for generating a vector
sequence x̃ ∈ `nT from available data z̃. For our main
theoretical results x̃ can be completely arbitrary, but it will be
shown below in Theorem 4 that best results are obtained if x̃
approximates a good state sequence for the model, hence we
refer to x̃ as the surrogate state sequence. Typical methods
for generating x̃ include subspace identification [29], using
a history of outputs as in nonlinear ARX [10], and – in the
context of model reduction – proper orthogonal decomposition
[2]. With some abuse of notation, from now on we consider
the available data set to include x̃, i.e. z̃ := {ũ, ỹ, x̃}.
Given a surrogate state sequence, we will frequently refer
to the following equation error functions:
(θ, z̃, t) = eθ (x̃(t + 1)) − fθ (x̃(t), ũ(t)),

η(θ, z̃, t) = ỹ(t) − gθ (x̃(t), ũ(t)),

(13)
(14)

although for brevity of notation we suppress the arguments
(θ, z̃, t) and simply denote dependence on time by subscripts:
t , ηt , and similarly for other signals e.g. xt = x(t).
Equation error t is the one-step prediction error in a coordinate system generated by the bijection e(·). A straightforward
way to search for a model is to minimize
PT the equation error in
a least-squares sense: JEE (θ, z̃) = t=0 |t |2 + |ηt |2 , which
is clearly convex if models are affinely parameterized. E.g.
[16] suggests identifying linear models in precisely this way.
However, while this is often effective there is in general no direct relationship between JEE and the (open-loop) simulation
error, which results from recursion as described in Section II.
Application of the Lagrangian relaxation (5) suggests the
following convex upper bound for simulation error, which
is valid for any choice of n-dimensional vector function
λ(x, z̃, t):
X

T
T
−1
X

2
0
sup
|g(xt , ũt ) − ỹt | +
λ e(xt+1 ) − f (xt , ũt ) .
n
x∈lT

t=0

t=0

Our first proposal for a convex upper bound on simulation
error comes from choosing λ(x, z̃, t) = −2(xt+1 − x̃t+1 ):
X
T
ˆ
JL (θ, z̃) := sup
|g(xt , ũt ) − ỹt |2
n
x∈lT

−

T
−1
X
t=0

t=0

0



2(xt+1 − x̃t+1 ) (e(xt+1 ) − f (xt , ũt )) ,

subject to x(0) = x̃(0).

(15)

With this upper bound and the model set from Definition 2,
we construct the following convex optimization problem.
Problem 1: For a model class (a, g, Θ) as in Definition 2
and data set z̃, solve JˆL? (z̃) = minθ∈Θ JˆL (θ, z̃).
A challenge in solving Problem 1 is that it involves computing the supremum of a function of a large number of variables,
the entire sequence x. An upper bound can be constructed
by interchanging the order of supremum and summation, and
introducing a positive-semidefinite storage function V (x, t):
J(θ, z̃) ≤

T
X
t=0

sup
xt ,xt+1

n
|g(xt , ũt ) − ỹt |2

+ λ(x, z̃, t)0 (e(xt+1 ) − f (xt , ũt ))
o
+ V (xt+1 , t + 1) − V (xt , t) .

(16)

Again, the quality of the bound will depend on the choice of
functions V and λ.
It will be shown below that a particular
PT choice of V and λ
leads to the upper bound JˆV (θ, z̃) := t=0 E(θ, z̃, t), where
E(θ, z̃,t) := sup {|f (x, ũt ) − f (x̃t , ũt ) − t |2(P −1 ) + |x − x̃t |2P
x∈Rn

− 2(x − x̃t )0 (e(x) − e(x̃t )) + |g(x, ũt ) − ỹt |2 } (17)
and t is given by (13). We define the following optimization
problem, which is convex in θ.
Problem 2: For a model class (a, g, Θ) as in Definition 2
and data set z̃, solve JˆV? (z̃) = minθ∈Θ JˆV (θ, z̃).
This upper bound is the robust identification error (RIE), first
introduced in [21]. The following definition characterizes data
sets for which “the true system is in the model class” [8].
Definition 4: Given a model class M̂ = (a, g, Θ) and a data
set z̃ we say that z̃ is in the range of M̂ if there exists a model
in M̂ satisfying (t) = η(t) = 0 for all t = 0, 1, ...T . That
is, simulating this model with the observed input and initial
conditions exactly reproduces the observed state and output.
We are now ready to state the main theoretical result
regarding Problems 1 and 2.
Theorem 4: Given a model class M̂ and a data set z̃, Let
?
(z̃) be the minimal simulation error over the model class
Jse
M̂, then for any data set z̃ and any model parameters θ in M̂,
Jse (θ, z̃) ≤ JˆL (θ, z̃) ≤ JˆV (θ, z̃) and minθ∈M̂ JˆV (θ, z̃) < ∞.
Furthermore, if the data set z̃ is in the range of M̂, then
?
Jse
(z̃) = JˆL? (z̃) = JˆV? (z̃) = 0.
Proof. The statement that Jse (θ, z̃) ≤ JˆL (θ, z̃) for all θ
follows from the basic properties of the Lagrangian relaxation.
To analyze JˆV (θ, z̃) we first show that it corresponds to
(16) with a particular choice of storage function and multiplier.
Indeed, let V (x, t) = 2(x − x̃t )0 (e(x) − e(x̃t )) − |x − x̃t |P and
λt = −2(xt+1 − x̃t+1 ), then (16) reduces to:
sup {|g(xt , ũt ) − ỹt |2 − |xt+1 − x̃t+1 |P

xt ,xt+1

+ 2(xt+1 − x̃t+1 )0 (f (xt , ũt ) − f (x̃t , ũt ) − t )

− 2(xt − x̃t )0 (e(xt ) − e(x̃t )) + |xt − x̃t |P }.

(18)

With this formulation, xt+1 only appears in quadratic or linear
terms, so the partial maximum can be computed exactly:
|f (xt , ũt ) − f (x̃t , ũt ) − t |2(P −1 ) = maxxt+1 {−|xt+1 − x̃t+1 |2P

+ 2(xt+1 − x̃t+1 )0 (f (xt , ũt ) − f (x̃t , ũt ) − t )} giving the
following bound in place of (18): supx(t) {|g(xt , ũt ) − ỹt |2 +
|f (xt , ũt ) − f (x̃t , ũt ) − t |2(P −1 ) − 2(xt − x̃t )0 (e(xt ) − e(x̃t )) +
|xt − x̃t |P } which is E(θ, z̃, t) in (17).
To prove the statement JˆL (θ, z̃) ≤ JˆV (θ, z̃), we first apply
(16) with the choice of V and λ given above, to give
X
T
sup
V (xt+1 , t + 1) − V (xt , t) + |g(xt , ũt ) − ỹt |2
x

t=0

+

λ0t (e(xt+1 )

− f (xt , ũt ))





≤ JˆV (θ, z̃).

(19)

All storage function terms cancel in the summation except
the first and the last. By construction2 , with x(0) = x̃(0)
we have V (x(0), 0) = 0, and since V is positive-definite
V (x(T ), 
T ) ≥ 0, the left-hand side of (19) is an upper bound


PT
2
0
for supx
|g(x
,
ũ
)−ỹ
|
+λ
(e(x
)−f
(x
,
ũ
))
t
t
t
t+1
t
t
t
t=0

which equals JˆL plus an additional multiplier term. This term
clearly does not decrease the supremum, since the choice
xt+1 = e−1 (f (xt , ũt )) would make it zero for any choice
of xt . This proves that JˆL (θ, z̃) ≤ JˆV (θ, z̃).
To prove that the minimal value of JˆV (θ, z̃) is finite, note
that one can choose a model (7) with f (x, u) = 0, g(x, u) =
0. In this case, the bound (17) reduces to E(θ, z̃, t) =
supx∈Rn {|ỹt |2 +|t |2(P −1 ) +|x−x̃t |2P −2(x−x̃t )0 (e(x)−e(x̃t ))}
and the stability constraint (10) ensures that |x − x̃t |2P − 2(x −
x̃t )0 (e(x) − e(x̃t )) ≤ −µ|x − x̃t |2 for all x, t. Since ỹt and t
are not functions of x this proves the supremum is finite.
Similarly, if z̃ is in the range of M̂, then there exists a
model θ̄ ∈ M̂ for which (t) = 0, η(t) = 0 for all t. Let
ē, f¯, ḡ be the model functions given by θ̄, then (17) is:
2

E(θ̄, z̃) = sup {|ḡ(x, ũt ) − ḡ(x̃t , ũt )| + |x −
x∈Rn

+ |f¯(x, ũt ) −

f¯(x̃t , ũt )|2(P −1 )

x̃t |2P

0

− 2(x − x̃t ) (ē(x) − ē(x̃t ))}

and by (10) this can be bounded by E(θ̄, z̃, t) ≤
supx∈Rn {−µ|x − x̃t |2 } = 0. Therefore JˆV (z̃, θ̄) = 0. Now,
since 0 ≤ Jse (θ, z̃) ≤ JˆL (θ, z̃) ≤ JˆV (θ, z̃) for all θ,
this model must give the minimal value of JˆV (z̃, θ̄), and
furthermore J(z̃, θ̄) = 0 so a perfect model is recovered. Note
that we do not require existence of a unique perfect model.
V. C ONTRACTING M ODELS AND L INEARIZED
S IMULATION E RROR
The upper bounds in the previous section, while convex,
may be expensive to compute when the dependence of e, f,
and g on x is nonlinear. Contraction analysis [30] provides an
alternative set of tools for studying model behaviour, based
on infinitesimal displacements of the system rather than pairs
of solutions. The principal advantage in the context of this
paper is that it results in significantly easier optimization
problems involving suprema of quadratic functions, which can
be computed explicitly or represented using LMIs.
2 Alternatively, when working without a known initial condition, one can
remove the constraint x(0) = x̃(0) and add −V (x0 , x̃0 ) = −2(x0 −
x̃0 )0 (e(x0 ) − e(x̃0 )) + |x0 − x̃0 |P to terms inside the brackets in (15).

A. Contraction Condition for Incremental `2 Stability
In this section we derive contraction conditions for incremental `2 stability. Contraction analysis is based on the study
of an extended system consisting of (8), (9) as well as the
differential dynamics
E(x(t + 1))∆(t + 1) = F (x(t), u(t))∆(t),

(20)

∆y (t) = G(x(t), u(t))∆(t),

(21)

i.e. the dynamics (8), (9) linearized with respect to initial
∂
∂
e(x), F (x, u) = ∂x
f (x, u), and
conditions, where E(x) = ∂x
∂
G(x, u) = ∂x g(x, u).
For our purposes, a differential storage function can be any
positive-semidefinite function of x and ∆, i.e. V : Rn ×Rn →
R such that V (x, ∆) ≥ 0 for all x and ∆. If a differential
storage function can be found that verifies the “differential”
dissipation inequality
V (x(t + 1), ∆(t + 1)) − V (x(t), ∆(t)) ≤
− |G(x(t), u(t))∆t |2 − µ|∆(t)|2

(22)

for solutions
is an upper
PTof (8), (20), then V (x(0),P∆(0))
T
bound for t=0 |G(x(t), u(t))∆t |2 + µ t=0 |∆(t)|2 for all
T , for all solutions of the extended system (8), (9), (20), (21).
We propose to use the Riemannian metric V (x, ∆) =
|E(x)∆|2(P −1 ) with P = P 0 > 0 as a differential storage
function, with which (22) becomes
|F (x, u)∆|2(P −1 ) − |E(x)∆|2(P −1 ) + |G(x, u)∆|2 ≤ −µ|∆|2 .
(23)
Just as (11) is implied by (10), (23) is implied by
F (x, u)0 P −1 F (x, u) + P
−E(x) − E(x)0 + G(x, u)0 G(x, u) ≤ −µI.

(24)

Furthermore, unlike (23), condition (24) is jointly convex in
e, f, g and P , i.e. convex in the parameter vector θ. We now
give an alternative model class to M̂ in Section III.
Definition 5: For a given choice of state dimension and
basis functions ei , fi , gi , we define Θ0 as the set of θ ∈ Rq
for which the model (8), (9) satisfies the inequality (24) for
some P = P 0 > 0. Additionally, we define the model class
0
M̂0 := (e−1
θ (fθ (x, u)), gθ (x, u), Θ ).
The main theoretical result for this model class is the
following, which is analogous to Theorems 1, 2, and 3.
Theorem 5: All models in M̂0 are well-posed and incrementally `2 stable. Addtionally, the set M̂0 contains all models
of the form (1), (2) for which (22) can be verified with
a constant quadratic metric V (x, ∆) = ∆0 M ∆ for some
M = M 0 > 0 and a(x, u) ∈ F.
Proof. To prove well-posedness, consider any pair x1 , x2 ∈
Rn and define xρ = (1 − ρ)x1 + ρx2 with ρ ∈ [0, 1]. Every
term on the left-hand side of (24) is positive semidefinite
except −E(x) − E(x)0 , soR E(x) + E(x)0 ≥ µI. Now,
1
(x2 −x1 )0 (e(x2 )−e(x1 )) = 0 (x2 −x1 )0 E(xρ )(x2 −x1 )dρ ≥
µ
2
2 |x2 − x1 | , so well-posedness is guaranteed by the same
argument as in the proof of Theorem 1.
To prove incremental stability, consider any pair of initial
conditions x1 (0) and x2 (0), and consider the family of solutions to (8) xρ , parameterized by ρ ∈ [0, 1], having initial

conditions xρ (0) = (1 − ρ)x1 (0) + ρx2 (0) and any fixed input
sequence ū(t). Note that for t > 0, xρ (t) does not necessarily
linearly interpolate x1 (t) and x2 (t).
∂
xρ (t). For any T > 0 we have
Now, define ∆ρ (t) := ∂ρ
2
PT
PT R 1
2
≤
t=0 |y2 (t) − y1 (t)| =
t=0 0 G(xρ (t))∆ρ (t)dρ
R 1 PT
2
t=0 |G(xρ (t))∆ρ (t)| dρ, by the Cauchy-Schwarz in0
equality and interchanging order of summation and integration.
Summing
(23)
repeatedly
with
x
=
xρ (t), u
=
ū(t) and R ∆
=
∆ρ (t) gives
PT
1
2
2
≤
|E(xρ (0))∆ρ (0)|(P −1 ) −
t=0 |y2 (t) − y1 (t)|
0
2
|E(xρ (T + 1))∆ρ (T + 1)|(P −1 ) dρ
R1
2
≤
|E(xρ (0))(x2 (0) − x1 (0))|(P −1 ) dρ, as ∆ρ (0) =
0
x2 (0) − x1 (0). As this bound holds uniformly for all T we
arrive at the desired
PTresult. Exactly the same reasoning can
be used to bound t=0 |x2 (t) − x1 (t)|2 .
The proof that all models with a constant quadratic metric
are in M̂0 is similar to the proof of Theorem 3, and follows
directly from the choice e(x) = M x, f (x, u) = M a(x, u).

~ = ~(θ, z̃) with the stacked vecrepresentation H(θ, z̃)∆
0 ˜
~
˜
˜ )0 ]0 and ~(θ, z̃) :=
tors ∆ := [∆(1) , ∆(2)0 , ..., ∆(T
0
0
0 0
[(0) , (1) , ...(T ) ] and H(θ, z̃) is a sparse matrix that is
affine in θ and is straightforward to construct from (28). The
analogue of (15) can then be written as

Condition (24) is, after taking the Schur complement of
F (x, u)0 P −1 F (x, u), a point-wise linear matrix inequality, i.e.
for every pair (x, u) there is an LMI. For polynomial models,
this can be enforced using a sum-of-squares constraint [31].

and the following convex cost function: JˆV0 (θ, z̃) =
PN
0
t=0 E (z̃, t) , which we call the local robust identification
error (local RIE).
Theorem 6: For all P = P 0 > 0 and signal data z̃, for every
model satisfying (24), J 0 (θ, z̃) ≤ JˆL0 (θ, z̃) ≤ JˆV0 (θ, z̃) < ∞.
Furthermore, if the data set z̃ is in the range of M̂, then the
optimal values satisfy J 0,? (z̃) = JˆL0,? (z̃) = JˆV0,? (z̃) = 0.
Proof. First we show that the bound JˆV0 is finite. Each term
E 0 (z̃) is the supremum of quadratic functional of ∆, so each
supremum is finite if and only if the corresponding matrix

B. Linearized Simulation Error and its Convex Bounds
We now introduce an approximation to simulation error that
is easier to optimize. Consider a “perturbed” version of the
model equations, motivated by (13), (14):
e(xρ (t + 1))

=

f (xρ (t), u(t)) + ρ(t),

(25)

yρ (t)

=

g(xρ (t), u(t)) + ρη(t),

(26)

where ρ ∈ [0, 1] and (t), η(t) are the equation errors.
By construction we have y0 (t) = y(t), the simulated
output, and y1 (t) = ỹ(t), the measured output. We define the
linearized simulation error as
N
X
1
|ỹ(t) − yρ (t)|2
ρ→1 (1 − ρ)2
t=0

J 0 (θ, z̃) := lim

(27)

to quantify local sensitivity of model equations with respect to equation errors. Applying L’Hospital’s rule twice
to (27) we arrive at the alternative expression J 0 (θ, z̃) =
PN
∂xρ (t)
2
˜
˜
t=0 |G(x̃(t), ũ(t))∆(t) + η(t)| , where ∆(t) =
∂ρ
˜
satisfies ∆(0)
= 0 and

˜ + 1) = F (x̃(t), ũ(t))∆(t)
˜
E(x̃(t + 1))∆(t
+ (t).

ρ=1

(28)

Roughly speaking, when the simulation error is small, i.e.
the model fit is good, the simulation error and the linearized
˜
simulation error will be very close, and ∆(t)
≈ x̃(t) − x(t).
Note that if the model equations are affine in x, then the simulation error and the linearized simulation error are identical.
A family of convex upper bounds for the linearized simulation error can be constructed in a similar way to those in
Section IV, so we skip some details in the construction.
The analogue of Problem 1 can be given in a simple form by noting that (28) can be written in “lifted”

~ ∆
~ + ~η (z̃)|2 − 2∆
~ 0 (H(θ, z̃)∆
~ − ~(z̃))}
JˆL0 (θ, z̃) := sup {|G(z̃)
~ nT
∆∈l
T

~
where G(z̃)
is a block-diagonal matrix with elements G(z̃(t))
and ~η is a stacked vector of output equation errors η(t). Note
that both H and ~ are affine in θ, so the functional JˆL0 is convex
in θ. Note also that the supremum is taken of a quadratic
~ which will be shown below to be concave, and
function of ∆,
can therefore be computed explicitly.
We can also construct an analogue of Problem 2 for linearized simulation error. Define the pointwise functions
E 0 (z̃, t) := sup {|F (x̃, ũ)∆ + |2(P −1 ) + |∆|2P
∆∈Rn

− 2∆0 E(x̃)∆ + |G(x̃, ũ)∆ + η|2 },

(29)

F (x̃(t), ũ(t))0 P −1 F (x̃(t), ũ(t)) + P
− E(x̃(t)) − E(x̃(t))0 + G(x̃(t), ũ(t))0 G(x̃(t), ũ(t))
is negative semidefinite. In fact, the contraction condition (23)
implies that the above matrix is is uniformly negative definite,
which implies a finite supremum achieved by a unique ∆.
The remainder of the proof follows from the same arguments as the proof of Theorem 4, but with the differential
storage function V (x̃, ∆) = ∆0 (E(x̃) + E(x̃)0 − P )∆ and by
˜ t for x̃t − xt , E(x̃t )∆
˜ t for e(x̃t ) − e(xt ) and
substituting ∆
˜
F (x̃t )∆t for f (x̃t , ut ) − f (xt , ut ).
VI. E XAMPLES AND D ISCUSSION
In this section we demonstrate the performance of the
proposed methods on two real-world applications: reducedorder dynamic inversion of an op-amp and identification of
a pneumatic actuator from experiments. Further examples
illustrating the method can be found in [24].
We judge the performance of a model over T time
samples
simulation error: Jperf =
qP using percent normalized
qP
PT
T
T
1
2
2
100
t=0 |ỹ(t) − y(t)| /
t=0 |ỹ(t) − T +1
s=0 ỹ(s)|
where y(t) is the result of simulating the model from
observed initial conditions and with the known input. When
minimizing the equation error, to ensure the dynamics are
well-posed, we require E(x) + E(x)0 − 2µI to be positive
semidefinite via a sums-of-squares constraint [31]. Software

TABLE II
C OMPARISON OF PERFORMANCE (JPERF ) FOR MODELS OF THE PNEUMATIC
PISTON ON VALIDATION DATA .
Model Class
Linear
Polynomial (degx (e)=1, degx (f )=1)
Polynomial (degx (e)=3, degx (f )=1)
Polynomial (degx (e)=3, degx (f )=3)

Equation Error
29.31
14.34
81.26
Divergent

Local RIE
24.62
22.07
14.84
6.64

Figure 2: Schematic of a two-stage op-amp.

(a) Simplified schematic of opamp.

(b) Opamp configuration.

Table 1: Design specifications

Fig. 1. Simplified internal schematic of the opamp and feedback configuration
openfor
loopdata
gain generation. Parasitics not
300 shown.
used
unity-gain bandwidth (ωu )
160MHz
phase margin
60◦
TABLE
I
CMRR
100
slew
rate
100MV/s
P ERFORMANCE (JPERF ) FOR O P -A MP √DYNAMIC
input-referred spot noise @ 1MHz (vnoise) 64nV/ Hz

INVERSION USING LOCAL

RIE WITH MODELS OF VARYING ORDER n AND DEGREE d.

n\d
1
2
3

1
86.47
23.82
11.80

3
94.162
21.90
8.31

5
94.22
21.69
7.99

packages such as the Yalmip toolbox [32] can be used to
transform polynomial optimization problems such as those
in this paper to semidefinite programs. The examples in this
circuitto
optimization
example .
paper Figure
took3: Aup
ten minutes
to solve using the free SDP
solver Sedumi [33], though preliminary work on specialised
biasing and performance
constraints, as illustrated
in the figure.
Designers canspeed-ups
also add constraints
algorithms
indicates
substantial
arebased
possible [34].

n their own design experiences. After solving this optimization problem, an optimal solution out of this

design space is obtained.
In this design example, with the specifications in table 1 and transistor modeled with 90nm predictive

A. Low-Order Nonlinear Dynamic Inversion of an Op-Amp
Our first example
is causal dynamic inversion of an op2
amp using data from a detailed high-dimensional simulation
model. Figure 1 presents the schematic of the op-amp at the
transistor level, and the device in the closed loop simulation
used for identification. For the task of dynamic inversion, we
seek to recover vin from vout . To accomplish this we identify
a dynamical system with vout as its input and vin as its output.
The identification data-set consists of two trials. Each trial is
a simulation in response to vin being the sum of five sinusoids.
The frequencies were chosen randomly within the bandwidth
of the amplifier (60 MHz). The trials have different biases and
amplitudes, and excite saturation in the amplifier.
Model complexity is varied in terms of state and input
dimension, and polynomial degree. Our input and surrogate state signals are in “autoregressive moving-average”
form, parameterized
ỹ(t) =
 by a positive integer n:
0
v
(t)
.
.
.
v
(t
+
1
−
n)
vin (t), ũ(t) =
, x̃(t) =
out
out

0
vin (t) . . . vin (t + 1 − n) . We search over models with
f being affine in the input and e, f of varying degree d in
the state x(t), with g fixed g(x(t), u(t)) = x1 (t). For all such
models, minimizing equation error JEE resulted in an unstable
model. Table I presents the results using RIE for varying model
complexities. The minimizers of the local RIE dramatically
improve in performance as the model complexity increases.
B. Identification of a Pneumatic Piston from Experiment
Our second example is the identification of a double-acting
pneumatic cylinder under conditions simulating heavy load.

Each chamber is separated from a control valve by approximately 25 cm of tubing. A compressor supplies high
pressure (between 0.6 and 0.8 MPa) which is then regulated to
a 0.5 MPa supply pressure to the valves, which have internal
pressure control mechanisms. Both ends of the piston are fixed
to a bracket. The dynamics of the system include changes in
supply pressure due to flow, filling the chambers, the internal
dynamics of the valves and elastic bending of the rail to which
the piston and bracket are fixed. The pressure is measured at
the inlet of both chambers (denoted p1 (t) and p2 (t) resp.). A
load cell between the bracket and piston measures the force
exerted on the bracket (fo (t)). We denote the commanded
pressure signals as v1 (t) and v2 (t) for valve 1 and 2 resp.
Finally, the supply pressure is measured near the inlet of the
valves (ps (t)). The input to the system is a scalar signal fd (t)
which represents the desired force in Newtons.
The experiment consisted of the application of seven chirp
signals of varying amplitudes being applied to the system.
These chirps signals swept linearly in instantaneous frequency
from approximately 3 to 40 Hz and peak desired force in
the range 20 to 100 N. The identification procedures were
trained on data sets with 20, 60, 100 and 130 N being the peak
absolute desired force. These models were then validated on
trials with 40, 80 and 120 N being the peak forces desired. The
experiments were performed at a sample-rate of 1 kHz. The
desired and realized force were taken to be the input and output
signals and the internal pressures and output force were chosen
0
as the surrogate states: ũ(t) = fd (t) fd (t − 1) , ỹ(t) =

0
fo (t), x̃(t) = p1 (t) p2 (t) ps (t) fo (t) . We fit several
discrete time models of the form (8), (9). We examine polynomial models of increasing complexity, beginning with a linear
model then input affine models with various degrees in x. In all
cases the output map g(·, ·) is fixed to g(x(t), u(t)) = x4 (t).
The computational complexity of minimizing the RIE scales
with the number of data-points for nonlinear fits. As a result,
we sub-select 10,000 data-points for optimizing the local
RIE. To ensure a fair comparison in these cases, we fit an
equation error model on both the same 10,000 data-points and
the whole data-set and report the better of the two scores.
Table II compares the results, in order of increasing model
complexity. Figure 2 compares the observed ỹ(t) to the open
loop simulation of the best performing model from the local
RIE. Note that the two plots are responses to differently scaled
versions of the same input, so nonlinearity is clearly apparent.
We observe that the local RIE under-performs equation
error for the first polynomial model structure. However, with
increasing model complexity the equation error minimization
at first improves, then rapidly degrades in performance and
becomes unstable. On the other hand, the performance of

150

Local RIE Model Performance for a Pneumatic Cylinder
Data
RIE Fit

Force (N)

100
50
0
−50
−100
−150

2

4

6

8

10

12

14

16

18

20

Time (s)
150

Data
RIE Fit

Force (N)

100
50
0
−50
−100
−150

2

4

6

8

10

12

14

16

18

20

Time (s)
(a) Data and Simulation.
Fig. 2. Model performance using the local RIE: simulation (solid-red) compared against observed output (gray dots) for two different input amplitudes
from nearly identical initial conditions. The response does not simply scale,
indicating nonlinearity.

the local RIE minimization steadily improves with increasing
model complexity. These observations are consistent with the
op-amp dynamic inversion.
One interpretation of this performance is that the stability
constraints and the “robust” fidelity bounds act as a form
of regularization, and mitigate the usual harmful effects of
over-fitting and reduce the need for structure selection or
“pruning” of regressors in autoregressive nonlinear models
[10]. Indeed, for linear systems we have observed that these
constraints confer a bias towards “more stable” models [22],
[34]. Quantifying this effect will be a topic of future research.
VII. ACKNOWLEDGMENTS
We would like to thank Ryuma Niiyama for assistance
with the pneumatic cylinder experiments, and Zhipeng Li, Yan
Li and Vladimir Marko Stojanovic for their collaboration in
providing the RF amplifier data set.
R EFERENCES
[1] A. C. Antoulas, Approximation of large-scale dynamical systems. Society for Industrial and Applied Mathematics, 2005, vol. 6.
[2] P. Holmes, J. Lumley, and G. Berkooz, Turbulence, Coherent Structures,
Dynamical Systems and Symmetry. Cambridge University Press, 1998.
[3] S. Lall, J. E. Marsden, and S. Glavaski, “A subspace approach to
balanced truncation for model reduction of nonlinear control systems,”
International Journal of Robust and Nonlinear Control, vol. 12, no. 6,
pp. 519–535, 2002.
[4] D. J. Lucia, P. S. Beran, and W. A. Silva, “Reduced-order modeling: new
approaches for computational physics,” Progress in Aerospace Sciences,
vol. 40, no. 1, pp. 51–117, 2004.
[5] 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.
[6] L. Ljung, System Identification: Theory for the User, 3rd ed. Englewood
Cliffs, New Jersey, USA: Prentice Hall, 1999.

[7] R. K. Pearson, “Nonlinear empirical modeling techniques,” Computers
& chemical engineering, vol. 30, no. 10, pp. 1514–1528, 2006.
[8] L. Ljung, “Perspectives on system identification,” Annual Reviews in
Control, vol. 34, no. 1, pp. 1 – 12, 2010.
[9] T. Schön, A. Wills, and B. Ninness, “System identification of nonlinear
state-space models,” Automatica, vol. 47, no. 1, pp. 39–49, 2010.
[10] S. Billings and S. Fakhouri, “Identification of systems containing linear
dynamic and static nonlinear elements,” Automatica, vol. 18, no. 1, pp.
15–26, 1982.
[11] S. Billings, Nonlinear System Identification: NARMAX Methods in the
Time, Frequency, and Spatio-Temporal Domains. Wiley, 2013.
[12] F. J. Doyle, R. K. Pearson, and B. A. Ogunnaike, Identification and
control using Volterra models. Springer Verlag, 2002.
[13] P. A. Regalia and P. Stoica, “Stability of multivariable least-squares
models,” IEEE Signal Processing Letters, vol. 2, no. 10, pp. 195–196,
1995.
[14] J. M. Maciejowski, “Guaranteed stability with subspace methods,”
Systems & Control Letters, vol. 26, no. 2, pp. 153–156, 1995.
[15] T. Van Gestel, J. A. Suykens, P. Van Dooren, and B. De Moor,
“Identification of stable models in subspace identification by using
regularization,” IEEE Transactions on Automatic Control, vol. 46, no. 9,
pp. 1416–1420, 2001.
[16] S. L. Lacy and D. S. Bernstein, “Subspace identification with guaranteed stability using constrained optimization,” IEEE Transactions on
Automatic Control, vol. 48, no. 7, pp. 1259–1263, 2003.
[17] B. Besselink, N. van de Wouw, and H. Nijmeijer, “Model reduction
for nonlinear systems with incremental gain or passivity properties,”
Automatica, vol. 49, no. 4, pp. 861–872, 2013.
[18] T. Söderström, “On the uniqueness of maximum likelihood identification,” Automatica, vol. 11, no. 2, pp. 193–197, 1975.
[19] A. Megretski, “Convex optimization in robust identification of nonlinear
feedback,” in Proceedings of the 47th IEEE Conference on Decision and
Control, Cancun, Mexico, Dec 9-11 2008, pp. 1370–1374.
[20] B. Bond, Z. Mahmood, Y. Li, R. Sredojevic, A. Megretski, V. Stojanovic,
Y. Avniel, and L. Daniel, “Compact modeling of nonlinear analog
circuits using system identification via semidefinite programming and
incremental stability certification,” IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems, vol. 29, no. 8, pp.
1149 –1162, aug. 2010.
[21] 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 IEEE Conference on Decision and Control
(CDC), 2010, pp. 7232–7237.
[22] I. R. Manchester, M. M. Tobenkin, and A. Megretski, “Stable nonlinear
system identification: Convexity, model class, and consistency,” in IFAC
Symposium on System Identification, vol. 16, no. 1, 2012, pp. 328–333.
[23] M. M. Tobenkin, I. R. Manchester, and A. Megretski, “Stable nonlinear
identification from noisy repeated experiments via convex optimization,”
in American Control Conference (ACC), 2013, pp. 3936–3941.
[24] M. M. Tobenkin, “Robustness analysis for identification and control
of nonlinear systems,” Ph.D. dissertation, Massachusetts Institute of
Technology, 2014.
[25] I. Pólik and T. Terlaky, “A survey of the S-lemma,” SIAM review, vol. 49,
no. 3, pp. 371–418, 2007.
[26] C. Lemaréchal, “Lagrangian relaxation,” in Computational Combinatorial Optimization. Springer, 2001, pp. 112–156.
[27] S. Boyd and L. Vandenberghe, Convex optimization.
Cambridge
university press, 2009.
[28] A. Poznyak, Advanced Mathematical Tools for Control Engineers:
Volume 1: Deterministic Systems. Elsevier, 2010.
[29] P. Van Overschee and B. De Moor, Subspace identification for linear
systems: theory, implementation, applications. Kluwer academic publishers, 1996.
[30] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for non-linear
systems,” Automatica, vol. 34, no. 6, pp. 683–696, June 1998.
[31] P. A. Parrilo, “Semidefinite programming relaxations for semialgebraic
problems,” Mathematical programming, vol. 96, pp. 293–320, 2003.
[32] J. Lofberg, “Yalmip: A toolbox for modeling and optimization in
matlab,” in Computer Aided Control Systems Design, 2004 IEEE International Symposium on. IEEE, 2004, pp. 284–289.
[33] J. F. Sturm, “Using SeDuMi 1.02, a matlab toolbox for optimization
over symmetric cones,” Optimization methods and software, vol. 11, no.
1-4, pp. 625–653, 1999.
[34] 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.

