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

TL;DR

This paper introduces convex optimization methods for nonlinear system identification and model reduction, improving stability and long-term prediction accuracy by leveraging Lagrangian relaxation, dissipation inequalities, and semidefinite programming.

Contribution

It presents novel convex optimization techniques for nonlinear system modeling, addressing non-convexity and local minima issues in traditional approaches.

Findings

Effective model order reduction for electronic circuits.

Successful pneumatic actuator identification from experimental data.

Enhanced stability and long-term prediction accuracy.

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 non-convex 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.

Tables

Table 1. TABLE I: Performance ( J perf subscript 𝐽 perf J_{\textrm{perf}} ) for Op-Amp dynamic inversion using local RIE with models of varying order n 𝑛 n and degree d 𝑑 d .

$n\setminus d$	1	3	5
1	86.47	94.162	94.22
2	23.82	21.90	21.69
3	11.80	8.31	7.99

Table 2. TABLE II: Comparison of performance ( J perf subscript 𝐽 perf J_{\textrm{perf}} ) for models of the pneumatic piston on validation data.

Model Class	Equation Error	Local RIE
Linear	29.31	24.62
Polynomial (degx($e$)=1, degx($f$)=1)	14.34	22.07
Polynomial (degx($e$)=3, degx($f$)=1)	81.26	14.84
Polynomial (degx($e$)=3, degx($f$)=3)	Divergent	6.64

Equations

$$x(t+1)$$

$$y(t)$$

$$J_{se}(\theta ,x,\tilde{z}):=\sum _{t=0}^T|y(t)-\tilde{y}(t){|}^2$$

$$y(t)=g(a(a(...a(x(0),\tilde{u}(0),\theta ),...,\tilde{u}(t-1),\theta ),\tilde{u}(t),\theta ),$$

$$J^{\star }=\underset{x,\theta }{\min }J(\theta ,x),\text{ subject to }h(\theta ,x)=0$$

$${\hat{J}}_{\lambda }^{\star }=\underset{\theta }{\min }{\hat{J}}_{\lambda }(\theta ),{\hat{J}}_{\lambda }(\theta ):=\underset{x}{\sup }\{J(\theta ,x)+\lambda (x{)}^{\prime }h(\theta ,x)\}$$

$$-c^{\prime }{P}^{-1}c\le b^{\prime }Pb-2b^{\prime }c$$

$$e_{\theta }(x)=\sum _{i=0}^q{\theta }_i{e}_i(x),f_{\theta }(x,u)=\sum _{i=0}^q{\theta }_i{f}_i(x,u),$$

$$g_{\theta }(x,u)=\sum _{i=0}^q{\theta }_i{g}_i(x,u),$$

$$e_{\theta }(x(t+1))$$

$$y(t)$$

$$|f(x_1,u)-f(x_2,u){|}_{(P^{-1})}^2-2(x_1-x_2{)}^{\prime }(e(x_1)-e(x_2))$$

$$+|x_1-x_2{|}_{P+\mu I}^2+|g(x_1,u)-g(x_2,u){|}^2\le 0$$

$$-|e(x_1)-e(x_2){|}_{(P^{-1})}^2+|f(x_1,u)-f(x_2,u){|}_{(P^{-1})}^2$$

$$+\mu |x_1-x_2{|}^2+|g(x_1,u)-g(x_2,u){|}^2$$

$$-|e(x_1)-e(x_2){|}_{(P^{-1})}^2\le |x_1-x_2{|}_P^2-2(x_1-x_2{)}^{\prime }(e(x_1)-e(x_2)).$$

$$|a(x_1,u)-a(x_2,u){|}_M^2-|x_1-x_2{|}_M^2\le -|g(x_1,u)-g(x_2,u){|}^2,$$

$$ϵ(\theta ,\tilde{z},t)$$

$$\eta (\theta ,\tilde{z},t)$$

$$\displaystyle\sup_{x\in l_{T}^{n}}\bigg{\{}$$

$${\hat{J}}_L(\theta ,\tilde{z}):=$$

$$\displaystyle-\sum_{t=0}^{T-1}2(x_{t+1}-\tilde{x}_{t+1})^{\prime}(e(x_{t+1})-f(x_{t},\tilde{u}_{t}))\bigg{\}},$$

$$\text{ subject to }x(0)=\tilde{x}(0).$$

$$J(\theta ,\tilde{z})\le$$

$$+\lambda (x,\tilde{z},t{)}^{\prime }(e(x_{t+1})-f(x_t,{\tilde{u}}_t))$$

$$\displaystyle+V(x_{t+1},t+1)-V(x_{t},t)\Big{\}}.$$

$$\mathcal{E}(\theta ,\tilde{z},$$

$$-2(x-{\tilde{x}}_t{)}^{\prime }(e(x)-e({\tilde{x}}_t))+|g(x,{\tilde{u}}_t)-{\tilde{y}}_t{|}^2\}$$

$$\underset{x_t,x_{t+1}}{\sup }\{|g(x_t,{\tilde{u}}_t)-{\tilde{y}}_t{|}^2-|x_{t+1}-{\tilde{x}}_{t+1}{|}_P$$

$$+2(x_{t+1}-{\tilde{x}}_{t+1}{)}^{\prime }(f(x_t,{\tilde{u}}_t)-f({\tilde{x}}_t,{\tilde{u}}_t)-{ϵ}_t)$$

$$-2(x_t-{\tilde{x}}_t{)}^{\prime }(e(x_t)-e({\tilde{x}}_t))+|x_t-{\tilde{x}}_t{|}_P\}.$$

$$\displaystyle\sup_{x}\bigg{\{}\sum_{t=0}^{T}$$

$$\displaystyle+\lambda_{t}^{\prime}(e(x_{t+1})-f(x_{t},\tilde{u}_{t}))\big{)}\bigg{\}}\leq\hat{J}_{V}(\theta,\tilde{z}).$$

$$\mathcal{E}(\bar{\theta },\tilde{z})=\underset{x\in {\mathbb{R}}^n}{\sup }\{|\bar{g}(x,{\tilde{u}}_t)-\bar{g}({\tilde{x}}_t,{\tilde{u}}_t){|}^2+|x-{\tilde{x}}_t{|}_P^2$$

$$+|\bar{f}(x,{\tilde{u}}_t)-\bar{f}({\tilde{x}}_t,{\tilde{u}}_t){|}_{(P^{-1})}^2-2(x-{\tilde{x}}_t{)}^{\prime }(\bar{e}(x)-\bar{e}({\tilde{x}}_t))\}$$

$$E(x(t+1))\Delta (t+1)=F(x(t),u(t))\Delta (t),$$

$${\Delta }_y(t)=G(x(t),u(t))\Delta (t),$$

$$V(x(t+1),\Delta (t+1))-V(x(t),\Delta (t))\le$$

Full text

I Introduction

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 model111We 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]., 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 (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-Impulse-Response (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 Problem Setup & Preliminaries

Spaces of vector signals, i.e. mappings $\{0,1,...,T\}\rightarrow\mathbb{R}^{n}$ are denoted $l_{T}^{n}$. For a symmetric matrix $A=A^{\prime}$, $A\geq 0$ and $A>0$ mean that $A$ is positive semidefinite and positive-definite, respectively. Inequalities $A\geq B,A>B$ mean $A-B\geq 0,A-B>0$, respectively. For symmetric positive-semidefinite $Q$ the notation $|v|_{Q}^{2}$ is shorthand for $v^{\prime}Qv$. Throughout the paper, signals with a tilde, e.g. $\tilde{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: $\tilde{u}=\{\tilde{u}(0),\tilde{u}(1),...,\tilde{u}(T)\}\in l_{T}^{m}$, $\tilde{y}=\{\tilde{y}(0),\tilde{y}(1),...,\tilde{y}(T)\}\in l_{T}^{p}$. We use the notation $\tilde{z}:=\{\tilde{u},\tilde{y}\}$ to represent the “available data”.

In this paper, a model class is a triple $\mathcal{M}=(a,g,\Theta)$ where $a:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{q}\rightarrow\mathbb{R}^{n}$ and $g:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{q}\rightarrow\mathbb{R}^{p}$ are continuously differentiable functions, and $\Theta\subset\mathbb{R}^{q}$ is a set of $q$-dimensional parameter vectors. A model class $(a,g,\Theta)$ is associated with state-space dynamical systems of the form

[TABLE]

where $x\in\mathbb{R}^{n},u\in\mathbb{R}^{m},y\in\mathbb{R}^{p}$, $\theta\in\Theta$. 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 $\ell^{2}$ stable if for every pair of initial conditions $x_{1}(0),x_{2}(0)$ and every input sequence $u(t),t=1,2,...$, the states $x_{1}(t),x_{2}(t)$ and outputs $y_{1}(t)=g(x_{1}(t),u(t))$ and $y_{2}(t)=g(x_{2}(t),u(t))$ satisfy $\sum_{t=0}^{\infty}|x_{1}(t)-x_{2}(t)|^{2}<\infty,\ \ \sum_{t=0}^{\infty}|y_{1}(t)-y_{2}(t)|^{2}<\infty.$

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:

[TABLE]

which is to be minimized over parameters $\theta$ and states $x$ for given data $\tilde{z}$, subject to dynamic constraints (1), (2).

The major difficulties we seek to address in this paper are:

1. (a)

Obtaining a useful set of stable models. For all non-trivial parameterizations of functions $a(x(t),u(t),\theta)$, an explicit and tractable characterization of the set of all $\theta$ defining a stable model is unknown. 2. (b)

Efficient global minimization of $J_{se}$. Even if the cost function $J(\theta,x,\tilde{z})$ is convex in $\theta$, the dynamic constraints (1) generally define a nonconvex set: the simulation output is generated by the recursion

[TABLE]

so even if $a$ and $g$ are linear in the parameter $\theta$, $J_{se}$ 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 $\hat{\mathcal{M}}=(a,g,\hat{\Theta})$ and optimization problems of the form: $\min_{\theta\in\hat{\Theta}}\hat{J}(\theta,\tilde{z}),$ with the following properties:

1. (a)

$\hat{\mathcal{M}}$ is a finite-dimensional convex parametrization of models satisfying the following inclusions: $LTI\subseteq\hat{\mathcal{M}}\subseteq\mathcal{M}_{\ell 2}$ where $LTI$ refers to the set of all stable linear time invariant models of appropriate dimension, and $\mathcal{M}_{\ell 2}$ is the set of all (possibly nonlinear) models that are incrementally $\ell^{2}$ stable. 2. (b)

$\hat{J}$ is a convex function of $\theta$ such that: (i) $\hat{J}(\theta,\tilde{z})\geq J_{se}(\theta,x,\tilde{z})$ for any feasible $x$; (ii) $\hat{J}(\theta,\tilde{z})$ always has a finite minimum over $\theta$; (iii) $\hat{J}(\theta,\tilde{z})=0$ when $\tilde{z}$ are generated by simulating the model represented by $\theta$. 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 S-Procedure, 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

[TABLE]

where $h(\theta,x)$ is a vector of constraint functions. In the context of simulation-error minimization, $J(\theta,x)$ represents the measure of model fit, e.g. (3), and $h(\theta,x)$ represents the dynamical constraints, e.g. (1).

Suppose that $J$ is a convex function of $\theta$ and $h$ is an affine function of $\theta$, but both are nonconvex functions of $\theta$ and $x$ jointly. In our context, Lagrangian relaxation refers to the approximation of the above problem by the following:

[TABLE]

where $\lambda$ is some vector function analogous to a Lagrange multiplier in nonlinear programming. The function $\hat{J}_{\lambda}(\theta)$ is clearly convex, since it is the supremum of a family of convex functions [27]. It is also an upper bound for $J(\theta,x)$ for any feasible $x$: if $h(\theta,x)=0$ in (5), then $\hat{J}_{\lambda}(\theta)=J(\theta,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: $h_{1}(\theta,x)=0\Leftrightarrow h_{2}(\theta,x)=0$, but which differ for other values of $\theta,x$, then the resulting values of $\hat{J}_{\lambda}(\theta)$ 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)=c^{\prime}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)=-c^{\prime}P^{-1}c$ with $P>0$ obey the following upper bound:

[TABLE]

where the right-hand side is a convex function of $c,P$ for any fixed $b$. Inequality (6) follows directly from the expansion $b^{\prime}Pb-2b^{\prime}c+c^{\prime}P^{-1}c=|c-Pb|^{2}_{P^{-1}}\geq 0.$ From this expansion it is also clear that the bound (6) is tight if $c=Pb$.

III Convex Parameterization of Stable Models

The models we propose are constructed from three linearly parametrized vector functions:

[TABLE]

where $\theta\in\mathbb{R}^{q}$ and each basis function $e_{i}:\ \mathbb{R}^{n}\to\mathbb{R}^{n},f_{i}:\ \mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n},g_{i}:\ \mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{p}$ is continuously differentiable in its arguments. These functions are associated with a parametrization of implicit state space dynamic models

[TABLE]

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,\theta)=e_{\theta}^{-1}(f_{\theta}(x,u))$ and $g(x,u,\theta)=g_{\theta}(x,u)$, as long as it can be guaranteed that $e_{\theta}(x)$ is a bijection for all $\theta$.

We denote by $\mathcal{E},\mathcal{F},\mathcal{G}$ the span of the functions (7) over $\theta\in\mathbb{R}^{q}$. We assume that $\mathcal{E},\mathcal{F},\mathcal{G}$ contain all linear functions $x\mapsto Ex,(x,u)\mapsto Ax+Bu,(x,u)\mapsto Cx+Du$, respectively, where $A,B,C,D,E$ are matrices of appropriate dimensions, and that $\mathcal{F}$ is closed under multiplication by a non-singular $n\times n$ matrix, i.e. if $f(\cdot,\cdot)\in\mathcal{F}$ then $Mf(\cdot,\cdot)\in\ F$ for all $M\in\mathbb{R}^{n\times n}$ with $\det(M)\neq 0$. Simple choices for $\mathcal{E},\mathcal{F},\mathcal{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 $\theta$ 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^{\prime}>0$:

[TABLE]

for all $x_{1},x_{2}\in\mathbb{R}^{n},u\in\mathbb{R}^{m}$ for some fixed $\mu>0$. Note that this constraint is convex with respect to $e,f,g,P$ and hence convex with respect to a parameter vector $\theta$ when $e,f,g,P$ are affinely parametrized.

Definition 2

For a given choice of state dimension and basis functions $e_{i},f_{i},g_{i}$, let $\Theta$ be the set of $\theta\in\mathbb{R}^{q}$ for which the model (8), (9) satisfies the inequality (10) for some $P=P^{\prime}>0$. Additionally, let $\hat{\mathcal{M}}$ be the model class $(e_{\theta}^{-1}(f_{\theta}(x,u)),g_{\theta}(x,u),\Theta)$.

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 $x_{1},x_{2},u$, then $e$ is a bijection $\mathbb{R}^{n}\leftrightarrow\mathbb{R}^{n}$.

*Proof. * We need to show that for every $z\in\mathbb{R}^{n}$, the equation $e(x)=z$ has a unique solution. Firstly, note that (10) implies $2(x_{1}-x_{2})^{\prime}(e(x_{1})-e(x_{2}))\geq\mu|x_{1}-x_{2}|^{2}.$ It is known (see, e.g., [28, Th. 18.15]) that if $T:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a nonlinear operator satisfying $(x_{1}-x_{2})^{\prime}(T(x_{1})-T(x_{2}))\geq c|x_{1}-x_{2}|^{2}$ for all $x_{1},x_{2}\in\mathbb{R}^{n}$ with $c>0$, then the nonlinear equation $T(x)=0$ has a unique solution. Take $T(x)=e(x)-z$ and $c=\mu/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 $\ell^{2}$ stability.

*Proof. * Consider the incremental Lyapunov function candidate $V(x_{1},x_{2})=|e(x_{1})-e(x_{2})|_{{(P^{-1})}}^{2}$ with $P=P^{\prime}>0$. The inequality

[TABLE]

for all $x_{1},x_{2}\in\mathbb{R}^{n},u\in\mathbb{R}^{m}$ is a sufficient condition for global incremental $\ell^{2}$ stability. Now apply the bound (6), with $c=e(x_{1})-e(x_{2})$ and $b=x_{1}-x_{2}$, giving

[TABLE]

and substitute the right-hand-side of this inequality for $-|e(x_{1})-e(x_{2})|^{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:\mathbb{R}^{n}\times\mathbb{R}^{m}\mapsto\mathbb{R}^{n}$ and $g:\mathbb{R}^{n}\times\mathbb{R}^{m}\mapsto\mathbb{R}^{p}$ is called quadratically incrementally $\ell^{2}$ stable, or more briefly: quadratically stable, if there exists an $M\in\mathbb{R}^{n\times n}$ with $M=M^{\prime}>0$ such that the condition

[TABLE]

holds for all $x_{1},x_{2}\in\mathbb{R}^{n}$ and $u\in\mathbb{R}^{m}$. I.e. incremental stability can be verified by a quadratic Lyapunov function.

Theorem 3

The model set $\hat{\mathcal{M}}$ includes all quadratically incrementally $\ell^{2}$ stable models of the form (1), (2) for which $a(\cdot,\cdot)\in\mathcal{F}$ and $g(\cdot,\cdot)\in\mathcal{G}$.

*Proof. * Suppose (12) holds for some $M=M^{\prime}>0$. Choose $e(x)=Mx,f(x,u)=Ma(x,u)$, and $P=M$, i.e. $c=Pb$ in (6), and leave $g(x,u)$ as is. Clearly this model is equivalent to (1), (2), and by assumption that $a(\cdot,\cdot)\in\mathcal{F}$ and that $\mathcal{F}$ is closed under multiplication by $M$, $f(\cdot,\cdot)\in\mathcal{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 Convex Upper Bounds for Simulation Error

Let us suppose there is method for generating a vector sequence $\tilde{x}\in\ell_{T}^{n}$ from available data $\tilde{z}$. For our main theoretical results $\tilde{x}$ can be completely arbitrary, but it will be shown below in Theorem 4 that best results are obtained if $\tilde{x}$ approximates a good state sequence for the model, hence we refer to $\tilde{x}$ as the surrogate state sequence. Typical methods for generating $\tilde{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 $\tilde{x}$, i.e. $\tilde{z}:=\{\tilde{u},\tilde{y},\tilde{x}\}$.

Given a surrogate state sequence, we will frequently refer to the following equation error functions:

[TABLE]

although for brevity of notation we suppress the arguments $(\theta,\tilde{z},t)$ and simply denote dependence on time by subscripts: $\epsilon_{t},\eta_{t}$, and similarly for other signals e.g. $x_{t}=x(t)$.

Equation error $\epsilon_{t}$ is the one-step prediction error in a coordinate system generated by the bijection $e(\cdot)$. A straightforward way to search for a model is to minimize the equation error in a least-squares sense: $J_{EE}(\theta,\tilde{z})=\sum_{t=0}^{T}|\epsilon_{t}|^{2}+|\eta_{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 $J_{EE}$ 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 $\lambda(x,\tilde{z},t)$:

[TABLE]

Our first proposal for a convex upper bound on simulation error comes from choosing $\lambda(x,\tilde{z},t)=-2(x_{t+1}-\tilde{x}_{t+1})$:

[TABLE]

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,\Theta)$ as in Definition 2 and data set $\tilde{z}$, solve $\hat{J}_{L}^{\star}(\tilde{z})=\min_{\theta\in\Theta}\hat{J}_{L}(\theta,\tilde{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)$:

[TABLE]

Again, the quality of the bound will depend on the choice of functions $V$ and $\lambda$.

It will be shown below that a particular choice of $V$ and $\lambda$ leads to the upper bound $\hat{J}_{V}(\theta,\tilde{z}):=\sum_{t=0}^{T}{\cal E}(\theta,\tilde{z},t),$ where

[TABLE]

and $\epsilon_{t}$ is given by (13). We define the following optimization problem, which is convex in $\theta$.

Problem 2

For a model class $(a,g,\Theta)$ as in Definition 2 and data set $\tilde{z}$, solve $\hat{J}_{V}^{\star}(\tilde{z})=\min_{\theta\in\Theta}\hat{J}_{V}(\theta,\tilde{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 $\hat{\mathcal{M}}=(a,g,\Theta)$ and a data set $\tilde{z}$ we say that $\tilde{z}$ is in the range of $\hat{\mathcal{M}}$ if there exists a model in $\hat{\mathcal{M}}$ satisfying $\epsilon(t)=\eta(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 $\hat{\mathcal{M}}$ and a data set $\tilde{z}$, Let $J_{se}^{\star}(\tilde{z})$ be the minimal simulation error over the model class $\hat{\mathcal{M}}$, then for any data set $\tilde{z}$ and any model parameters $\theta$ in $\hat{\mathcal{M}}$, $J_{se}(\theta,\tilde{z})\leq\hat{J}_{L}(\theta,\tilde{z})\leq\hat{J}_{V}(\theta,\tilde{z})$ and $\min_{\theta\in\hat{\mathcal{M}}}\hat{J}_{V}(\theta,\tilde{z})<\infty$. Furthermore, if the data set $\tilde{z}$ is in the range of $\hat{\mathcal{M}}$, then $J_{se}^{\star}(\tilde{z})=\hat{J}_{L}^{\star}(\tilde{z})=\hat{J}_{V}^{\star}(\tilde{z})=0.$

*Proof. * The statement that $J_{se}(\theta,\tilde{z})\leq\hat{J}_{L}(\theta,\tilde{z})$ for all $\theta$ follows from the basic properties of the Lagrangian relaxation.

To analyze $\hat{J}_{V}(\theta,\tilde{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-\tilde{x}_{t})^{\prime}(e(x)-e(\tilde{x}_{t}))-|x-\tilde{x}_{t}|_{P}$ and $\lambda_{t}=-2(x_{t+1}-\tilde{x}_{t+1})$, then (16) reduces to:

[TABLE]

With this formulation, $x_{t+1}$ only appears in quadratic or linear terms, so the partial maximum can be computed exactly: $|f(x_{t},\tilde{u}_{t})-f(\tilde{x}_{t},\tilde{u}_{t})-\epsilon_{t}|_{(P^{-1})}^{2}=\max_{x_{t+1}}\{-|x_{t+1}-\tilde{x}_{t+1}|^{2}_{P}\\ +2(x_{t+1}-\tilde{x}_{t+1})^{\prime}(f(x_{t},\tilde{u}_{t})-f(\tilde{x}_{t},\tilde{u}_{t})-\epsilon_{t})\}$ giving the following bound in place of (18): $\sup_{x(t)}\{|g(x_{t},\tilde{u}_{t})-\tilde{y}_{t}|^{2}+|f(x_{t},\tilde{u}_{t})-f(\tilde{x}_{t},\tilde{u}_{t})-\epsilon_{t}|_{(P^{-1})}^{2}-2(x_{t}-\tilde{x}_{t})^{\prime}(e(x_{t})-e(\tilde{x}_{t}))+|x_{t}-\tilde{x}_{t}|_{P}\}$ which is ${\cal E}(\theta,\tilde{z},t)$ in (IV).

To prove the statement $\hat{J}_{L}(\theta,\tilde{z})\leq\hat{J}_{V}(\theta,\tilde{z})$, we first apply (16) with the choice of $V$ and $\lambda$ given above, to give

[TABLE]

All storage function terms cancel in the summation except the first and the last. By construction222Alternatively, when working without a known initial condition, one can remove the constraint $x(0)=\tilde{x}(0)$ and add $-V(x_{0},\tilde{x}_{0})=-2(x_{0}-\tilde{x}_{0})^{\prime}(e(x_{0})-e(\tilde{x}_{0}))+|x_{0}-\tilde{x}_{0}|_{P}$ to terms inside the brackets in (15)., with $x(0)=\tilde{x}(0)$ we have $V(x(0),0)=0$, and since $V$ is positive-definite $V(x(T),T)\geq 0$, the left-hand side of (19) is an upper bound for \sup_{x}\bigg{\{}\sum_{t=0}^{T}\big{(}|g(x_{t},\tilde{u}_{t})-\tilde{y}_{t}|^{2}+\lambda_{t}^{\prime}(e(x_{t+1})-f(x_{t},\tilde{u}_{t}))\big{)}\bigg{\}} which equals $\hat{J}_{L}$ plus an additional multiplier term. This term clearly does not decrease the supremum, since the choice $x_{t+1}=e^{-1}(f(x_{t},\tilde{u}_{t}))$ would make it zero for any choice of $x_{t}$. This proves that $\hat{J}_{L}(\theta,\tilde{z})\leq\hat{J}_{V}(\theta,\tilde{z})$.

To prove that the minimal value of $\hat{J}_{V}(\theta,\tilde{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 (IV) reduces to ${\cal E}(\theta,\tilde{z},t)=\sup_{x\in\mathbb{R}^{n}}\{|\tilde{y}_{t}|^{2}+|\epsilon_{t}|_{(P^{-1})}^{2}+|x-\tilde{x}_{t}|_{P}^{2}-2(x-\tilde{x}_{t})^{\prime}(e(x)-e(\tilde{x}_{t}))\}$ and the stability constraint (10) ensures that $|x-\tilde{x}_{t}|_{P}^{2}-2(x-\tilde{x}_{t})^{\prime}(e(x)-e(\tilde{x}_{t}))\leq-\mu|x-\tilde{x}_{t}|^{2}$ for all $x,t$. Since $\tilde{y}_{t}$ and $\epsilon_{t}$ are not functions of $x$ this proves the supremum is finite.

Similarly, if $\tilde{z}$ is in the range of $\hat{\mathcal{M}}$, then there exists a model $\bar{\theta}\in\hat{\mathcal{M}}$ for which $\epsilon(t)=0,\eta(t)=0$ for all $t$. Let $\bar{e},\bar{f},\bar{g}$ be the model functions given by $\bar{\theta}$, then (IV) is:

[TABLE]

and by (10) this can be bounded by ${\cal E}(\bar{\theta},\tilde{z},t)\leq\sup_{x\in\mathbb{R}^{n}}\{-\mu|x-\tilde{x}_{t}|^{2}\}=0.$ Therefore $\hat{J}_{V}(\tilde{z},\bar{\theta})=0$. Now, since $0\leq J_{se}(\theta,\tilde{z})\leq\hat{J}_{L}(\theta,\tilde{z})\leq\hat{J}_{V}(\theta,\tilde{z})$ for all $\theta$, this model must give the minimal value of $\hat{J}_{V}(\tilde{z},\bar{\theta})$, and furthermore $J(\tilde{z},\bar{\theta})=0$ so a perfect model is recovered. Note that we do not require existence of a unique perfect model.

V Contracting Models and Linearized Simulation Error

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.

V-A Contraction Condition for Incremental $\ell^{2}$ Stability

In this section we derive contraction conditions for incremental $\ell^{2}$ stability. Contraction analysis is based on the study of an extended system consisting of (8), (9) as well as the differential dynamics

[TABLE]

i.e. the dynamics (8), (9) linearized with respect to initial conditions, where $E(x)=\frac{\partial}{\partial x}e(x)$, $F(x,u)=\frac{\partial}{\partial x}f(x,u)$, and $G(x,u)=\frac{\partial}{\partial x}g(x,u)$.

For our purposes, a differential storage function can be any positive-semidefinite function of $x$ and $\Delta$, i.e. $V:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ such that $V(x,\Delta)\geq 0$ for all $x$ and $\Delta$. If a differential storage function can be found that verifies the “differential” dissipation inequality

[TABLE]

for solutions of (8), (20), then $V(x(0),\Delta(0))$ is an upper bound for $\sum_{t=0}^{T}|G(x(t),u(t))\Delta_{t}|^{2}+\mu\sum_{t=0}^{T}|\Delta(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,\Delta)=|E(x)\Delta|_{(P^{-1})}^{2}$ with $P=P^{\prime}>0$ as a differential storage function, with which (22) becomes

[TABLE]

Just as (11) is implied by (10), (23) is implied by

[TABLE]

Furthermore, unlike (23), condition (24) is jointly convex in $e,f,g$ and $P$, i.e. convex in the parameter vector $\theta$. We now give an alternative model class to $\hat{\mathcal{M}}$ in Section III.

Definition 5

For a given choice of state dimension and basis functions $e_{i},f_{i},g_{i}$, we define $\Theta^{0}$ as the set of $\theta\in\mathbb{R}^{q}$ for which the model (8), (9) satisfies the inequality (24) for some $P=P^{\prime}>0$. Additionally, we define the model class $\hat{\mathcal{M}}^{0}:=(e_{\theta}^{-1}(f_{\theta}(x,u)),g_{\theta}(x,u),\Theta^{0})$.

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 $\hat{\mathcal{M}}^{0}$ are well-posed and incrementally $\ell^{2}$ stable. Addtionally, the set $\hat{\mathcal{M}}^{0}$ contains all models of the form (1), (2) for which (22) can be verified with a constant quadratic metric $V(x,\Delta)=\Delta^{\prime}M\Delta$ for some $M=M^{\prime}>0$ and $a(x,u)\in\mathcal{F}$.

*Proof. * To prove well-posedness, consider any pair $x_{1},x_{2}\in\mathbb{R}^{n}$ and define $x_{\rho}=(1-\rho)x_{1}+\rho x_{2}$ with $\rho\in[0,1]$. Every term on the left-hand side of (24) is positive semidefinite except $-E(x)-E(x)^{\prime}$, so $E(x)+E(x)^{\prime}\geq\mu I$. Now, $(x_{2}-x_{1})^{\prime}(e(x_{2})-e(x_{1}))=\int_{0}^{1}(x_{2}-x_{1})^{\prime}E(x_{\rho})(x_{2}-x_{1})d\rho\geq\frac{\mu}{2}|x_{2}-x_{1}|^{2}$, 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 $x_{1}(0)$ and $x_{2}(0)$, and consider the family of solutions to (8) $x_{\rho}$, parameterized by $\rho\in[0,1]$, having initial conditions $x_{\rho}(0)=(1-\rho)x_{1}(0)+\rho x_{2}(0)$ and any fixed input sequence $\bar{u}(t)$. Note that for $t>0$, $x_{\rho}(t)$ does not necessarily linearly interpolate $x_{1}(t)$ and $x_{2}(t)$.

Now, define $\Delta_{\rho}(t):=\frac{\partial}{\partial\rho}x_{\rho}(t)$. For any $T>0$ we have

$\sum_{t=0}^{T}|y_{2}(t)-y_{1}(t)|^{2}=\>\sum_{t=0}^{T}\left|\int_{0}^{1}G(x_{\rho}(t))\Delta_{\rho}(t)d\rho\right|^{2}\leq\int_{0}^{1}\sum_{t=0}^{T}\left|G(x_{\rho}(t))\Delta_{\rho}(t)\right|^{2}d\rho$, by the Cauchy-Schwarz inequality and interchanging order of summation and integration.

Summing (23) repeatedly with $x=x_{\rho}(t)$, $u=\bar{u}(t)$ and $\Delta=\Delta_{\rho}(t)$ gives $\sum_{t=0}^{T}|y_{2}(t)-y_{1}(t)|^{2}\leq\int_{0}^{1}\left|E(x_{\rho}(0))\Delta_{\rho}(0)\right|_{(P^{-1})}^{2}-\left|E(x_{\rho}(T+1))\Delta_{\rho}(T+1)\right|_{(P^{-1})}^{2}\;d\rho\\ \leq\int_{0}^{1}\left|E(x_{\rho}(0))(x_{2}(0)-x_{1}(0))\right|_{(P^{-1})}^{2}d\rho,$ as $\Delta_{\rho}(0)=x_{2}(0)-x_{1}(0)$. As this bound holds uniformly for all $T$ we arrive at the desired result. Exactly the same reasoning can be used to bound $\sum_{t=0}^{T}|x_{2}(t)-x_{1}(t)|^{2}$.

The proof that all models with a constant quadratic metric are in $\hat{\mathcal{M}}^{0}$ is similar to the proof of Theorem 3, and follows directly from the choice $e(x)=Mx,f(x,u)=Ma(x,u)$.

Condition (24) is, after taking the Schur complement of $F(x,u)^{\prime}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].

V-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):

[TABLE]

where $\rho\in[0,1]$ and $\epsilon(t),\eta(t)$ are the equation errors.

By construction we have $y_{0}(t)=y(t)$, the simulated output, and $y_{1}(t)=\tilde{y}(t)$, the measured output. We define the linearized simulation error as

[TABLE]

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}(\theta,\tilde{z})=\sum_{t=0}^{N}|G(\tilde{x}(t),\tilde{u}(t))\tilde{\Delta}(t)+\eta(t)|^{2},$ where $\tilde{\Delta}(t)=\left.\frac{\partial x_{\rho}(t)}{\partial\rho}\right|_{\rho=1}$ satisfies $\tilde{\Delta}(0)=0$ and

[TABLE]

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 $\tilde{\Delta}(t)\approx\tilde{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” representation $H(\theta,\tilde{z})\vec{\Delta}=\vec{\epsilon}(\theta,\tilde{z})$ with the stacked vectors $\vec{\Delta}:=[\tilde{\Delta}(1)^{\prime},\tilde{\Delta}(2)^{\prime},...,\tilde{\Delta}(T)^{\prime}]^{\prime}$ and $\vec{\epsilon}(\theta,\tilde{z}):=[\epsilon(0)^{\prime},\epsilon(1)^{\prime},...\epsilon(T)^{\prime}]^{\prime}$ and $H(\theta,\tilde{z})$ is a sparse matrix that is affine in $\theta$ and is straightforward to construct from (28). The analogue of (15) can then be written as

[TABLE]

where $\vec{G}(\tilde{z})$ is a block-diagonal matrix with elements $G(\tilde{z}(t))$ and $\vec{\eta}$ is a stacked vector of output equation errors $\eta(t)$. Note that both $H$ and $\vec{\epsilon}$ are affine in $\theta$, so the functional $\hat{J}^{0}_{L}$ is convex in $\theta$. Note also that the supremum is taken of a quadratic function of $\vec{\Delta}$, which will be shown below to be concave, and can therefore be computed explicitly.

We can also construct an analogue of Problem 2 for linearized simulation error. Define the pointwise functions

[TABLE]

and the following convex cost function: $\hat{J}^{0}_{V}(\theta,\tilde{z})=\sum_{t=0}^{N}{\cal E}^{0}(\tilde{z},t)$ , which we call the local robust identification error (local RIE).

Theorem 6

For all $P=P^{\prime}>0$ and signal data $\tilde{z}$, for every model satisfying (24), $J^{0}(\theta,\tilde{z})\leq\hat{J}^{0}_{L}(\theta,\tilde{z})\leq\hat{J}^{0}_{V}(\theta,\tilde{z})<\infty.$ Furthermore, if the data set $\tilde{z}$ is in the range of $\hat{\mathcal{M}}$, then the optimal values satisfy $J^{0,\star}(\tilde{z})=\hat{J}_{L}^{0,\star}(\tilde{z})=\hat{J}_{V}^{0,\star}(\tilde{z})=0.$

*Proof. * First we show that the bound $\hat{J}^{0}_{V}$ is finite. Each term ${{\cal E}}^{0}(\tilde{z})$ is the supremum of quadratic functional of $\Delta$, so each supremum is finite if and only if the corresponding matrix

[TABLE]

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 $\Delta$.

The remainder of the proof follows from the same arguments as the proof of Theorem 4, but with the differential storage function $V(\tilde{x},\Delta)=\Delta^{\prime}(E(\tilde{x})+E(\tilde{x})^{\prime}-P)\Delta$ and by substituting $\tilde{\Delta}_{t}$ for $\tilde{x}_{t}-x_{t}$, $E(\tilde{x}_{t})\tilde{\Delta}_{t}$ for $e(\tilde{x}_{t})-e(x_{t})$ and $F(\tilde{x}_{t})\tilde{\Delta}_{t}$ for $f(\tilde{x}_{t},u_{t})-f(x_{t},u_{t})$.

VI Examples and Discussion

In this section we demonstrate the performance of the proposed methods on two real-world applications: reduced-order 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 using percent normalized simulation error: $J_{\textrm{perf}}=100{\sqrt{\sum_{t=0}^{T}|\tilde{y}(t)-y(t)|^{2}}}/{\sqrt{\sum_{t=0}^{T}|\tilde{y}(t)-\frac{1}{T+1}\sum_{s=0}^{T}\tilde{y}(s)|^{2}}}$ 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)^{\prime}-2\mu I$ to be positive semidefinite via a sums-of-squares constraint [31]. Software 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 paper took up to ten minutes to solve using the free SDP solver Sedumi [33], though preliminary work on specialised algorithms indicates substantial speed-ups are possible [34].

VI-A Low-Order Nonlinear Dynamic Inversion of an Op-Amp

Our first example is causal dynamic inversion of an op-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 $v_{\textrm{in}}$ from $v_{\textrm{out}}$. To accomplish this we identify a dynamical system with $v_{\textrm{out}}$ as its input and $v_{\textrm{in}}$ as its output.

The identification data-set consists of two trials. Each trial is a simulation in response to $v_{\textrm{in}}$ 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 by a positive integer $n$: $\tilde{y}(t)=v_{\textrm{in}}(t),\>\tilde{u}(t)=\begin{bmatrix}v_{\textrm{out}}(t)&\ldots&v_{\textrm{out}}(t+1-n)\end{bmatrix}^{\prime},\>\tilde{x}(t)=\begin{bmatrix}v_{\textrm{in}}(t)&\ldots&v_{\textrm{in}}(t+1-n)\end{bmatrix}^{\prime}.$ 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))=x_{1}(t)$. For all such models, minimizing equation error $J_{EE}$ 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.

VI-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 $p_{1}(t)$ and $p_{2}(t)$ resp.). A load cell between the bracket and piston measures the force exerted on the bracket ($f_{o}(t)$). We denote the commanded pressure signals as $v_{1}(t)$ and $v_{2}(t)$ for valve 1 and 2 resp. Finally, the supply pressure is measured near the inlet of the valves ($p_{s}(t)$). The input to the system is a scalar signal $f_{d}(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 as the surrogate states: $\tilde{u}(t)=\begin{bmatrix}f_{d}(t)&f_{d}(t-1)\end{bmatrix}^{\prime},\quad\tilde{y}(t)=f_{o}(t),\tilde{x}(t)=\begin{bmatrix}p_{1}(t)&p_{2}(t)&p_{s}(t)&f_{o}(t)\end{bmatrix}^{\prime}.$ 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(\cdot,\cdot)$ is fixed to $g(x(t),u(t))=x_{4}(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 $\tilde{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 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.