On the smoothness of nonlinear system identification

TL;DR

This paper analyzes the smoothness challenges in nonlinear system identification and proposes multiple shooting as a solution to improve optimization stability and feasibility in complex models.

Contribution

It introduces multiple shooting for nonlinear system identification, addressing smoothness issues and enabling estimation in chaotic or unstable systems.

Findings

Lipschitz constant and smoothness can grow exponentially with simulation length in non-contractive regions.

Multiple shooting improves optimization stability by splitting data into smaller subsets.

The method effectively estimates parameters of chaotic systems and neural networks.

Abstract

We shed new light on the \textit{smoothness} of optimization problems arising in prediction error parameter estimation of linear and nonlinear systems. We show that for regions of the parameter space where the model is not contractive, the Lipschitz constant and $\beta$-smoothness of the objective function might blow up exponentially with the simulation length, making it hard to numerically find minima within those regions or, even, to escape from them. In addition to providing theoretical understanding of this problem, this paper also proposes the use of multiple shooting as a viable solution. The proposed method minimizes the error between a prediction model and the observed values. Rather than running the prediction model over the entire dataset, multiple shooting splits the data into smaller subsets and runs the prediction model over each subset, making the simulation length a design parameter and making it possible to solve problems that would be infeasible using a standard approach. The equivalence to the original problem is obtained by including constraints in the optimization. The new method is illustrated by estimating the parameters of nonlinear systems with chaotic or unstable behavior, as well as neural networks. We also present a comparative analysis of the proposed method with multi-step-ahead prediction error minimization.

Figures

Tables

Table 1. Table 1: Computational cost for the logistic map estimation. The number of function evaluations and total running time until convergence. Minimum, maximum and median are given over 15 runs for the situations presented in Fig.  2 . *The number of iterations is limited to 1000 and the solver is interrupted when this number is reached.

	function evaluations			run time (s)
$\Delta m_{\max }$	min	median	max	min	median	max
$N$	1	15	23	0.01	0.2	1.1
$10$	43	1000*	1000*	0.7	24.7	27.3
$5$	29	115	645	1.0	2.9	26.9
$2$	21	50	65	1.7	2.9	3.8

Table 2. Table 2: Tradeoffs: multiple shooting vs MSA-PEM. N θ subscript 𝑁 𝜃 N_{\theta} is the number of parameters. N x subscript 𝑁 𝑥 N_{x} is the number of states that need to be propagated through the simulation ( N x = n y = 2 subscript 𝑁 𝑥 subscript 𝑛 𝑦 2 N_{x}=n_{y}=2 ). Δ ​ m max Δ subscript m \Delta{\rm m}_{\max} is the maximum simulation length used for multiple shooting and K 𝐾 K is the number of steps ahead predicted in the MSA-PEM.

	multiple shoot.	MSA-PEM
Problem dimension	$N_{\theta }+\frac{N}{\Delta m_{\max }}\cdot N_x$	$N_{\theta }$
# of constraints	$\frac{N}{\Delta m_{\max }}\cdot N_x-1$	$\varnothing$
Time-complexity (computing $V$)	$𝒪\left(N\left(n_y+n_u\right)\right)$	$𝒪\left(NK\left(n_y+n_u\right)\right)$
Asymp. statistical properties	Independent of $\Delta m_{\max }$	Approach expected properties as $K\to N$
# of independent subproblems	$\frac{N}{\Delta m_{\max }}$	$N$

Equations

$$V=\frac{1}{N}\sum _{k=1}^N\Vert \mathbf{y}[k]-\hat{\mathbf{y}}[k]{\Vert }^2.$$

$$\underline{\mathbf{u}}[k]$$

$$\underline{\mathbf{y}}[k-1]$$

$$\mathbf{y}[k]={\mathbf{f}}^{\ast }(\underline{\mathbf{y}}[k-1],\underline{\mathbf{u}}[k])+\mathbf{v}[k],$$

$$\hat{\mathbf{y}}[k]={\mathbf{f}}_{\mathbf{\theta }}(\underline{\mathbf{y}}[k-1],\underline{\mathbf{u}}[k]).$$

$$\bar{\mathbf{y}}[k]$$

$$\mathbf{y}[k]$$

$$\tilde{\mathbf{y}}[k]$$

$$\hat{\mathbf{y}}[k]$$

$$\mathbf{y}[k]={\mathbf{f}}^{\ast }(\underline{\mathbf{y}}[k-1],\underline{\mathbf{u}}[k],\mathbf{v}[k-1],\dots ,\mathbf{v}[k-n_v])+\mathbf{v}[k].$$

$$\tilde{\mathbf{v}}[k]=\mathbf{y}[k]-{\mathbf{f}}_{\mathbf{\theta }}(\underline{\mathbf{y}}[k-1],\underline{\mathbf{u}}[k],\tilde{\mathbf{v}}[k-1],\dots ,\tilde{\mathbf{v}}[k-n_v])$$

$$\hat{\mathbf{y}}[k]={\mathbf{f}}_{\mathbf{\theta }}(\underline{\mathbf{y}}[k-1],\underline{\mathbf{u}}[k],\tilde{\mathbf{v}}[k-1],\dots ,\tilde{\mathbf{v}}[k-n_v]).$$

$$\mathbf{x}[k]$$

$$\hat{\mathbf{y}}[k]$$

$$\underset{\mathbf{\theta },{\mathbf{x}}_0}{\min }V.$$

$$L_V=\{\begin{array}{ll}\mathcal{O}(L_h^{2N})& \text{if }{L}_h>1,\\ \mathcal{O}(N)& \text{if }{L}_h=1,\\ \mathcal{O}(1)& \text{if }{L}_h<1.\end{array}$$

$$L_V^{\prime }=\{\begin{array}{ll}\mathcal{O}(L_h^{3N})& \text{if }{L}_h>1,\\ \mathcal{O}\left(N^3\right)& \text{if }{L}_h=1,\\ \mathcal{O}\left(1\right)& \text{if }{L}_h<1.\end{array}$$

$${\mathbf{x}}^i[k]$$

$${\hat{\mathbf{y}}}^i[k]$$

$$V_i=\frac{1}{\Delta {\mathrm{m}}_i}\sum _{k={\mathrm{m}}_i+1}^{{\mathrm{m}}_{i+1}}\Vert \mathbf{y}[k]-{\hat{\mathbf{y}}}^i[k]{\Vert }^2,$$

$$V^M=\sum _{i=1}^M\frac{\Delta {\mathrm{m}}_i}{N}{V}_i.$$

$$\underset{\mathbf{\theta },{\mathbf{x}}_0^1,\cdots ,{\mathbf{x}}_0^M}{\min }$$

$${\mathbf{x}}^{i-1}[{\mathrm{m}}_i]={\mathbf{x}}_0^i,$$

$$\text{for }i=2,3,\cdots ,M.$$

$$L_{V^M}=\mathcal{O}(L_h^{2\Delta {\mathrm{m}}_{\max }});L_{V^M}^{\prime }=\mathcal{O}(L_h^{3\Delta {\mathrm{m}}_{\max }}).$$

$$y[k]=\theta y[k-1](1-y[k-1]),$$

$$\displaystyle\left\{\begin{array}[]{l}x_{1}[k+1]=x_{1}[k]+\delta x_{2}[k]\\ x_{2}[k+1]=-\delta\frac{g}{l}\sin x_{1}[k]+(1-\delta\frac{k_{a}}{m})x_{2}[k]+\delta\frac{1}{m}u[k]\end{array}\right.$$

$$y[k]=x_1[k]$$

$$\bar{y}[k]$$

$$y[k]$$

$$D_{\mathbf{\theta }}[k]$$

$$J_{\mathbf{\theta }}[k]$$

$$D_{{\mathbf{x}}_0}[k]$$

$$J_{{\mathbf{x}}_0}[k]$$

$$\nabla V=\frac{2}{N}\sum _{k=1}^{N}J[k{]}^T(\hat{\mathbf{y}}[k]-\mathbf{y}[k]),$$

$${\nabla }^{2}V=\frac{2}{N}\sum _{k=1}^N\left(J[k{]}^{T}J[k]+\mathbf{S}[k]\right).$$

Full text

, , , ,

1 Introduction

Prediction error methods [1] are a widespread class of methods for parameter estimation of dynamic models, which estimate the parameters by minimizing the error between predicted and measured trajectories. Many well-known estimation methods fit into this framework, such as minimizing the one-step-ahead prediction error or the free-run-simulation error. While the classical literature focuses primarily on the estimation of linear systems [1], the framework is general and enjoys appealing asymptotic properties for the general nonlinear setup [2].

Minimizing the one-step-ahead prediction usually yields an easier optimization problem, the minimization of the free-run simulation error or other recurrent structures, however, may produce more accurate models. These recurrent models often have a smaller generalization error [3], [4] and better capability when it comes to long-term prediction [5]. Minimizing recurrent structures is used, for instance, to improve the model structure selection of polynomial models [6], fine-tune parameters of nonlinear state-space [7] and block-oriented models [8].

It is common knowledge among practitioners that the optimization problem resulting from a recurrent model structure is harder to solve [9]. For linearly parametrized models and convex loss functions, minimization of one-step-ahead prediction error leads to a convex optimization problem; for recurrent model structures, the ensuing optimization is, in general, non-convex, complicating the search for global optima. Even during local optimization, recurrent model structures can lead to cost functions with poor smoothness properties, including many ‘jagged’ local minima, cf. Fig. 2(a) for an illustration. The current understanding of the relationship between model internal dynamics and the smoothness properties of the cost function is, however, imprecise and provides little insight into ways of circumventing the problem. Furthermore, the few studies that do investigate the objective function properties in this context are focused on linear systems, see e.g. [10].

The purpose of this paper is twofold. First, we aim to provide insight into the properties of the objective function arising in prediction error estimation problems in a general nonlinear setup. Specifically, we show how the smoothness of the objective function depends on two factors: the simulation length and the decay rate of the recurrent part of the prediction model. Second, we illustrate how this theoretical insight might be leveraged for the design and analysis of practical system identification methods.

The use of multiple shooting is analyzed in the context of prediction error minimization. This technique reformulates the optimization problem that arises from minimizing the difference between the output of a prediction model and the observed values. Rather than running the prediction model over the entire dataset, the multiple shooting formulation splits the dataset into smaller subsets and runs the prediction model over each subset. The equivalence with the original problem is obtained by including equality constraints in the optimization problem. This method results in a smoother objective function since it works with shorter simulations and prevents trajectories from diverging too much.

The multiple shooting formulation has reportedly provided improvements in the parameter estimation of ordinary differential equations [11], [12], [13] and in the solution of optimal control [14], [15], [16]. In the context of system identification, multiple shooting has been used for estimating polynomial nonlinear space-state models [17] and output error models [18] in settings where conventional methods fail to provide good solutions. Here, we extend this method to the entire class of prediction error methods. In addition, theoretical arguments are put forward to help understand why and when the proposed method is useful. We also present a comparative analysis with multi-step-ahead prediction error minimization [19], [20] [21], showing strengths and weaknesses of each method at a conceptual level and, also, with numerical examples.

2 Prediction error methods

While prediction error methods are widely known, they are often introduced from a linear perspective [1]. In the present section, this is accomplished in a nonlinear setting.

Consider the dataset $\mathcal{Z}^{N}=\{(\mathbf{u}[k],\mathbf{y}[k]),k=1,2,\cdots,N\}$ containing $N$ measured inputs and outputs of a dynamical system. Prediction error methods assume an internal model that delivers output predictions $\hat{\mathbf{y}}[k]$, $k=1,2,\cdots,N$. A cost function is defined as the distance between predictions and measured values:

[TABLE]

The prediction model is, usually, parametrized by a parameter vector $\boldsymbol{\theta}$ and, although such a dependence is not made explicit in the notation, $\hat{\mathbf{y}}[k]$ depends upon $\boldsymbol{\theta}$. An estimate $\hat{\boldsymbol{\theta}}$ of the parameter vector may be obtained by minimizing $V$.

Next we assume a dynamical, stochastic, discrete-time system as the data generating process. Let:

[TABLE]

where $n_{y}$, $n_{u}$ are the maximum input and output lags. The examples we give next will differ in how they account for the noise in the model. But, in a noiseless situation, the data generating process for all of them corresponds to a difference equation $\mathbf{y}[k]=\mathbf{f}^{*}(\underline{\mathbf{y}}[k-1],\underline{\mathbf{u}}[k])$. The parametrized function $\mathbf{f}_{\boldsymbol{\theta}}$ is defined and, if the noise assumption and model structure are correct, the minimization of the cost function $V$ would yield the estimate $\hat{\boldsymbol{\theta}}$ such that, as $N\rightarrow\infty$, then $\mathbf{f}_{\mathbf{\hat{\boldsymbol{\theta}}}}\rightarrow\mathbf{f}^{*}$ (or to a function with different structure but the same performance in predicting the training dataset). See Appendix D for a more precise description of the asymptotic properties of nonlinear prediction error methods.

Choosing the “true” model structure (one for which there exists a $\boldsymbol{\theta}^{*}$ such that $\mathbf{f}_{\boldsymbol{\theta}^{*}}=\mathbf{f}^{*}$) might be impossible in a practical application. Nevertheless the assumption is not so restrictive as it might appear, since there exist families of universal approximator functions (e.g. neural networks and polynomials) for which the distance ${\|\mathbf{f}_{\boldsymbol{\theta}}-\mathbf{f}^{*}\|}$ might be made arbitrarily small within a compact set.

2.1 Nonlinear ARX models

The nonlinear ARX (autoregressive with exogenous input) model encodes an output that is corrupted by white process noise. That is, it assumes the data were generated by the stochastic discrete-time system:

[TABLE]

where $\mathbf{v}[k]$ is a zero-mean white noise. This assumption yields (cf. Appendix D) the prediction model:

[TABLE]

The minimization of the cost function in Eq. (1) for this prediction model yields an estimator with the desired asymptotic properties.

2.2 Nonlinear OE models

The OE (output error) model encodes an output that is corrupted by white measurement noise. That is, it is assumed that the data were generated by:

[TABLE]

where, again, $\mathbf{v}[k]$ is zero-mean white noise and $\bar{\mathbf{y}}[k]$ represents the noiseless output. This assumption yields the following prediction model:

[TABLE]

Here $\tilde{\mathbf{y}}[k]$ represents an estimate of the noiseless output $\bar{\mathbf{y}}[k]$, which should approach the true value as ${\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}^{*}}$.

2.3 Nonlinear ARMAX models

The nonlinear ARMAX (autoregressive moving average with exogenous input) model encodes an output that is corrupted by additive zero-mean process noise. In this case, the propagation equation allows the noise term $\mathbf{v}[k]$ to be propagated by the dynamics:

[TABLE]

This assumption allows the model to account for some forms of colored process noise and results in the following prediction model:

[TABLE]

Here $\tilde{\mathbf{v}}[k]$ represents an estimate of the noise corrupting the system and will approach the true noise if the estimated parameter vector approaches the true parameter value $\boldsymbol{\theta}^{*}$.

2.4 General nonlinear state-space framework

From now on, we will focus on a state-space representation that is general enough to encompass the prediction model from the above examples (for an appropriate choice of the functions $\mathbf{h}$ and $\mathbf{g}$). For this representation, the predicted output is given by:

[TABLE]

where $\textbf{x}[k]$ denotes the internal state vector at instant $k$. For the ARX model the transition state would be empty $\mathbf{x}=\varnothing$; for the OE model ${\mathbf{x}[k]=\left(\bar{\mathbf{y}}[k-1],\cdots,\bar{\mathbf{y}}[k-n_{v}]\right)}$; and, for the ARMAX model $\mathbf{x}[k]=\left(\tilde{\mathbf{v}}[k-1],\cdots,\tilde{\mathbf{v}}[k-n_{y}]\right)$.

Here ${\underline{\mathbf{z}}[k]=(\underline{\mathbf{y}}[k-1],\underline{\mathbf{u}}[k])}$. Using both inputs $\underline{\mathbf{u}}[k]$ and autoregressive terms $\underline{\mathbf{y}}[k-1]$ is what allows this state-space representation to encompass ARX, ARMAX and OE models, and also other representations such as the polynomial greybox models proposed in [22].

2.5 Initial conditions

In order to guarantee the desirable asymptotic properties, the prediction model needs to be simulated starting with the appropriate initial conditions $\mathbf{x}_{0}$. Since the true initial condition $\mathbf{x}_{0}^{*}$ is unknown, there are two possible approaches when estimating the parameters.

The first approach is to fix $\mathbf{x}_{0}$, for some ${\mathbf{x}_{0}\approx\mathbf{x}_{0}^{*}}$, and minimize the cost function (1). This approach is based on the idea that, for an asymptotically stable system, the influence of the initial conditions on the output will decrease over time for many cases of interest [23] and, hence, even if $\mathbf{x}_{0}\not=\mathbf{x}_{0}^{*}$ we can still obtain a good estimate of the parameters. In this case the first samples may be discarded, to make sure the transient errors are not too large. For the ARMAX model, an appropriate choice of initial values would be ${\tilde{\mathbf{v}}[k]=\mathbf{0},~{}k=1,\cdots,n_{v}}$ and, for the OE model, ${\tilde{\mathbf{y}}[k]=\mathbf{y}[k],~{}k=1,\cdots,n_{y}}$.

The second approach consists in including $\mathbf{x}_{0}$ in the optimization problem, so it converges to $\mathbf{x}_{0}^{*}$ and improves the quality of the parameter estimates. The optimization problem to be solved in this case is to minimize $V$ with both $\boldsymbol{\theta}$ and $\mathbf{x}_{0}$ as optimization variables:

[TABLE]

3 Smoothness of prediction error methods

The theorem below relates the Lipschitz constant of $V$, and its gradient (i.e. $\beta$-smoothness), to the simulation length $N$. The Lipschitz constant of the cost function and the $\beta$-smoothness both play a crucial role in optimization [24]. Lower values imply that local (Taylor) expansions of the cost function are more predictive of the cost function, and that the optimization algorithm can still converge while taking larger steps. It also gives an upper bound on how distinct in performance two close local minima may be.

The first part of the theorem below can be seen as a formalization of the exploding gradient problem, often studied in the context of neural networks [25]. The second part provides information about the explosion of second-order derivatives and curvature and it is, to the best of our knowledge, novel. As a result of the analysis, it is found that not only walls (resulting from large first-order derivatives) might be formed in non-contractive regions of the parameter space, but also regions with exploding curvature with multiple close local minima (cf. Fig. 2(a)). A recurrent neural network-oriented perspective of these results is discussed in a concurrent work from our group [26].

Theorem 1**.**

Let $\textbf{h}(\mathbf{x},\underline{\mathbf{z}};\boldsymbol{\theta})$ and $\textbf{g}(\mathbf{x},\underline{\mathbf{z}};\boldsymbol{\theta})$ in (6) be Lipschitz in ${(\mathbf{x},\boldsymbol{\theta})}$ with constants $L_{h}$ and $L_{g}$ on a compact and convex set $\Omega=(\Omega_{\mathbf{x}},\Omega_{\underline{\mathbf{z}}},\Omega_{\boldsymbol{\theta}})$. With $\{\underline{\mathbf{z}}[k]\}_{k=1}^{N}\subseteq\Omega_{\underline{\mathbf{z}}}$ and $(\Omega_{\mathbf{x}},\Omega_{\boldsymbol{\theta}})\subseteq\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{\theta}}$. If there exists at least one choice of ${(\mathbf{x}_{0},\boldsymbol{\theta})}$ for which there is an invariant set contained in $\Omega$, then, for trajectories and parameters confined within $\Omega$:

The cost function $V$ defined in (1) is Lipschitz with constant:111Where $\mathcal{O}$ denotes the big O notation. It should be read as: $L(N)=\mathcal{O}(g(N))$ if and only if there exist positive integers $M$ and $N_{0}$ such that ${|L(N)|\leq\;Mg(N)\text{ for all }N\geq N_{0}}$.

[TABLE] 2. 2.

If the Jacobian matrices of h and g are also Lipschitz with respect to ${(\mathbf{x},\boldsymbol{\theta})}$ on $\Omega$, then the gradient of the cost function $\nabla V$ is also Lipschitz with constant:

[TABLE]

{pf}

See Appendix B.

For contractive models222We say that a dynamical system $\mathbf{x}[k+1]=\mathbf{h}(\mathbf{x}[k])$ is contractive if, for all $\mathbf{x}$ and $\mathbf{w}$, it satisfies $\|\mathbf{h}(\mathbf{x})-\mathbf{h}(\mathbf{w})\|<L\|\mathbf{x}-\mathbf{w}\|$, for $L<1$., under certain regularity conditions, we have $L_{h}<1$ and, according to the above theorem, both the Lipschitz constant and the $\beta$-smoothness of the cost function can be bounded by a constant that, asymptotically, does not depend on the simulation length. All contractive systems have a unique fixed point inside the contractive region, and all trajectories converge to such a fixed point [27, Theorem 9.23]. Systems with richer nonlinear dynamic behaviors, such as limit cycles and chaotic attractors, and also unstable systems, are non-contractive and will always have $L_{h}\geq 1$. The Lipschitz constants and $\beta$-smoothness for these systems may, according to Theorem 1, blow up exponentially (or polynomially for some limit cases) with the maximum simulation length.

Less formally, for models that have infinitely long dependencies (i.e. are non-contractive) the distance between trajectories of models that are close in the parameter space might become progressively larger along the simulation length because errors will accumulate. This might yield very intricate objective functions in some parts of the parameter space and make the optimization problem very dependent on the initial estimate.

4 Multiple shooting

The theoretical results from the previous section suggest that long simulation lengths might yield regions of the parameter space where the cost function is intricate, hence hard for the optimization algorithm to navigate.

In this section, we propose the application, in the context of prediction error methods, of a technique called multiple shooting for which the maximum simulation length is a design parameter. This enables solving problems that would be impossible or very hard to solve in the setting of Section 2, which will be named single shooting from now on.

4.1 Method formulation

For the multiple shooting formulation, rather than simulating the prediction model (6) through the entire dataset from a single initial condition vector $\mathbf{x}_{0}$, the data is split into $M$ intervals ${\{[{\rm m}_{i}+1,{\rm m}_{i+1}]\mid i=1,\cdots,M\}}$, $0={\rm m}_{1}<{\rm m}_{2}<\cdots<{\rm m}_{M}<{\rm m}_{M+1}=N$, each one with its own set of initial conditions $\mathbf{x}_{0}^{i}\in\mathbb{R}^{N_{x}}$. The $i$-th vector of initial conditions $\mathbf{x}_{0}^{i}$ is used to compute $\hat{\mathbf{y}}^{i}[k]$ over ${{\rm m}_{i}+1\leq k\leq{\rm m}_{i+1}}$:

[TABLE]

Since the length of the simulation is limited to the shorter interval $[{\rm m}_{i}+1,{\rm m}_{i+1}]$, the trajectory is less likely to strongly diverge and this typically helps the optimization procedure by making the objective function smoother.

Let ${\Delta{\rm m}_{i}={\rm m}_{i+1}-{\rm m}_{i}}$, we define:

[TABLE]

to be the cost function associated with the $i$-th interval, where the prediction $\hat{\mathbf{y}}^{i}[k]$ depends upon $\boldsymbol{\theta}$ and $\mathbf{x}_{0}^{i}$, according to (10). The multiple shooting formulation makes use of the following objective function:

[TABLE]

This objective function includes states $\mathbf{x}_{0}^{1},\cdots,\mathbf{x}_{0}^{M}$ as free variables in the optimization. Hence, rather than reinforcing the cohesion of the states $\mathbf{x}[k]$ by defining them through a recurrence relation that casts a dependency of $\mathbf{x}[k]$ all the way back to the initial condition $\mathbf{x}_{0}$, as in the single shooting formulation, the cohesion between subsequent states is achieved through optimization constraints, resulting in the following problem:

[TABLE]

The next theorem provides the equivalence between (7) and (12). The theorem and its corollary are a formalization of the intuition provided in Fig. 1 and they give further insight into how the constraints in the multiple shooting formulation are used to imitate a single simulation throughout the entire dataset.

Theorem 2**.**

If $\mathbf{x}^{i-1}[{\rm m}_{i}]=\mathbf{x}_{0}^{i},\text{for }i=2,3,\cdots,M$ and $\mathbf{x}_{0}^{1}=\mathbf{x}_{0}$, then $V=V^{M}$.

{pf}

Let us call $\mathbf{x}[k]$, $\hat{\mathbf{y}}[k]$ and $\mathbf{x}^{i}[k]$, $\hat{\mathbf{y}}^{i}[k]$ the states and predictions in, respectively, the single shooting simulation and in the $i$-th multiple shooting interval. For a fixed $i$, if $\mathbf{x}[{\rm m}_{i}]=\mathbf{x}_{0}^{i}$ then $\mathbf{x}[k]=\mathbf{x}^{i}[k]$ for all ${k\in[{\rm m}_{i}+1,{\rm m}_{i+1}]}$. Hence, inside this same interval, $\hat{\mathbf{y}}[k]=\hat{\mathbf{y}}^{i}[k]$. Applying this for every $i$ it follows from the respective definitions that: $V=\sum_{i=1}^{M}\tfrac{\Delta m_{i}}{N}~{}V_{i}=V^{M}$.

Corollary 3**.**

The pair $(\boldsymbol{\theta}^{*},\mathbf{x}_{0}^{*})$ is a global solution of (7) if and only if there exist $(\mathbf{x}_{0}^{2},\cdots,\mathbf{x}_{0}^{M})$ such that $(\boldsymbol{\theta}^{*},\mathbf{x}_{0}^{*},\mathbf{x}_{0}^{2},\cdots,\mathbf{x}_{0}^{M})$ is a global solution of the optimization problem (13).

Multiple shooting can be understood as a generalization of the single shooting case. That is because, if $M=1$ and $\Delta{\rm m}_{1}=N$, both methods result in the same optimization problem. Multiple shooting, however, might provide some advantages. The method is more amenable to parallelization, since each cost function $V_{i}$ and its associated derivatives can be computed independently and, possibly, in parallel. Also, long simulations usually yield larger numerical errors due to finite precision errors that accumulate along the simulation, hence multiple shooting shorter simulation intervals also attenuate this problem. Finally, multiple shooting cost function has better smoothness properties that will be investigated next.

4.2 Properties of the cost function

For the multiple shooting method, the Lipschitz constant of $V^{M}$ and of its gradient, $L_{V^{M}}$ and $L_{V^{M}}^{\prime}$ depend asymptotically on $\Delta\text{m}_{\max}=\max_{1\leq i\leq M}\Delta\text{m}_{i}$ rather than on $N$ (See Appendix C). For instance, if $L_{h}>1$:

[TABLE]

Since $\Delta\rm{m}_{\max}$ is a design parameter, it is possible to have some control over the Lipschitzness and $\beta$-smoothness of the objective function in the non-contractive region of the parameter space ($L_{h}\geq 1$).

5 Implementation and numerical examples

In this section, numerical examples are presented. The cost function smoothness is investigated through the lens of the theoretical results in Section 3 and it is shown how multiple-shooting might help mitigate some problems.

The equality-constrained problem (12), which arises from the multiple shooting formulation, is solved using an implementation of the sequential quadratic programming solver originally described in [28] available in the SciPy library [29]333 scipy.optimize.minimize(method=‘trust-constr’). The procedure used for computing the derivatives is explained in Appendix A. Additional numerical examples are provided in Appendix F.

5.1 OE model for a chaotic system

This example illustrates how multiple shooting makes prediction error methods more robust w.r.t. the choice of initial conditions for the optimization. A dataset with $N=200$ samples is generated using the logistic map [30]:

[TABLE]

with $\theta=3.78$. From the generated dataset we try to estimate an output error model with the same structure.

Fig. 2(a) illustrates the objective function for the single shooting case. The data generating system, Eq. (15), presents chaotic behavior for $\theta\in[3.57,4]$, which explains the very intricate objective function in this region. For chaotic systems, small variations in the parameters may cause large variations in the system trajectory and, consequently, abrupt changes in the free-run simulation error. This explains why the cost function used in estimating an OE model has many local minima for this problem.

The solutions found by the solver, for different initial guesses, are also displayed in Fig. 2(a). Notice that the solver fails to find the true solution because it always gets trapped in local stationary points near the initial guess. Even in this noiseless scenario, the estimation problem is very challenging due to the chaotic nature of the system. Hence, for a simulation that is sufficiently long, the trajectories will differ significantly even for small parameter variations.

Multiple shooting makes the problem easier by limiting the maximum simulation length. Fig. 2 (b), (c) and (d) display the objective function and the solutions found by the solver starting from different initial guesses. The estimation procedure becomes easier as $\Delta\rm{m}_{\max}$ is made smaller. For Fig. 2(c) the solver already converges to the true parameter for some initial guesses, but not for all of them. For Fig. 2(d), the solver converges to the true solution for all initial guesses that have been tested.

For the multiple shooting case, besides $\theta$, the initial conditions are also optimization parameters. To help with the visualization of this multidimensional problem, Fig. 2(b), (c) and (d) display the main curve corresponding to the objective function for the true initial conditions and faded lines corresponding to the objective function for perturbed initial conditions. Another consequence of the problem having more parameters than displayed in the figure is that the cost function found by the solver does not need to lie on any of the objective function curves displayed in the figure, since it may have a different set of initial conditions $\mathbf{x}_{0}^{i}$.

It is important to highlight that the mechanism used here is not to take the system outside of the chaotic regime, but rather avoid simulating the system for too long. By doing that, we avoid the major problem that arises in the identification of chaotic systems, i.e. the high sensitivity to parameters and initial conditions. This results in a better behaved objective function (cf. Fig. 2). The constraints allow the equivalence with the original prediction error problem (according to Theorem 2 and Corollary 3).

Table 1 gives the number of function evaluations and the running time for the four situations displayed in Fig. 2. The convergence happens within just a few iterations for $\Delta{\rm m}_{\max}=N$ (single shooting) because any initial point is probably very close to some optimal local solution. As we reduce $\Delta\rm{m}_{\max}$ the objective function becomes less intricate and this is reflected in the convergence of the solver. For $\Delta\rm{m}_{\max}=10$ the solver takes much longer to converge. We believe this happens because the local solution is not so close in the parameter space to the initial guess anymore. As $\Delta\rm{m}_{\max}$ is further decreased, however, the convergence becomes faster because it is dealing with, what is believed to be, a smoother problem that can be handled more accurately by low order approximations.

5.2 Pendulum and inverted pendulum

Consider the following discrete-time nonlinear system:

[TABLE]

which corresponds to a pendulum model, discretized using the Euler approximation ${\dot{x}(t)\approx\frac{x((k+1)\delta)-x(k\delta)}{\delta}}$, where $g$ is the gravity acceleration, $m$ is the mass connected to the extremity of the pendulum, $l$ is the length of the (massless) rod connecting the mass to the pivot point, and $k_{a}$ is the linear friction constant. It has two states: the angle of the mass ($x_{1}$) and the angular velocity ($x_{2}$). The input $u[k]$ is the force applied to the mass.

This system has multiple equilibrium points, namely, $(x_{1},x_{2})=(\pm\pi i,0)$ for $i=0,1,2,3,\cdots$. The equilibrium points at $(x_{1},x_{2})=(\pm 2\pi i,0)$ are stable and the equilibrium points at $(x_{1},x_{2})=(\pi\pm 2\pi i,0)$ are unstable. For this system, with $g=9.8$, $l=0.3$, $m=3$, $k_{a}=2$ and $\delta=0.01$, we define three different datasets: (a) A dataset for which small inputs are applied to the system, that stays under the influence of the stable point $(x_{1},x_{2})=(0,0)$ and $y[k]$ stays, approximately, inside the range $\left[-\frac{\pi}{2},+\frac{\pi}{2}\right]$; (b) A dataset for which the system is maintained close to the unstable point $(x_{1},x_{2})=(\pi,0)$ by a linear controller; and, (c) A dataset for which the input is large enough to drive the pendulum to full rotations around its center. The output corresponding to those three situations are displayed in Fig. 3.

Fixing $m=3$ and $\delta=0.01$ parameters $\frac{g}{l}$ and $k_{a}$ of an output error model with the structure presented in (19) were estimated from the data. A visualization of the cost function is presented in Fig. 4 together with numerical solutions found by means of the single shooting and multiple shooting formulation starting from different initial conditions.

For dataset (a), the single shooting formulation is able to recover the true parameters from data for most of the initial conditions. Some exceptions occur when initialized far away from the correct initial conditions. For datasets (b) and (c), for which the system needs, respectively, to operate close to the unstable dynamics or to account for the existence of multiple fixed points, the cost function is highly intricate and full of local minima. In this case, the optimization algorithm, even when initialized close to the local solution, fails to converge to reasonable solutions. This result is consistent with Theorem 1 and how the smoothness of the objective function degenerates (exponentially) on sets of the parameter space for which the prediction model is non-contractive, such as the trajectories close to the unstable fixed point of the system (19). The use of multiple shooting yields an objective function that looks similar to a paraboloid in the region of interest for the three cases and the optimization procedure converges to the true parameter regardless of the initialization.

For nonlinear ARX, ARMAX or OE models the states are directly measured (possibly with some noise contamination) and the initialization of the optimization parameters $\mathbf{x}_{0}^{i}$ follows naturally from the considerations in Section 2.5. Here, however, the state variable $x_{2}$ is not measured. This variable can be interpreted as the derivative of $x_{1}$, so a finite difference approximation was used in the initialization of the intermediary initial conditions $\mathbf{x}_{0}^{i}$. Although the high-pass behavior of the derivative amplifies the noise, this choice was still better than a completely arbitrary one.

6 Comparison with multi-step-ahead prediction error minimization

6.1 Multi-step-ahead prediction error minimization

Multiple-shooting is presented here as a possible way of limiting the simulation interval $\Delta\rm{m}_{\max}$. A method that appears in the system identification literature that also has a similar effect is the multi-step-ahead prediction error minimization (MSA-PEM) [19], [20] [21]. The approach fits well into the moving horizon framework [20] and is popular for system identification in model predictive control application. The method has been studied primarily in a linear model setting, nevertheless it can be extended to a nonlinear setting [21].

In the MSA-PEM estimation, the simulation is truncated to a fixed number $K$ of steps backwards. That is, for each $k=1,\cdots,N$, we define an auxiliary variable $\tilde{\mathbf{x}}_{k}[i]$ and an initial condition $\tilde{\mathbf{x}}_{0,k}$. Starting from $\tilde{\mathbf{x}}_{k}[k-K]=\tilde{\mathbf{x}}_{0,k}$ the system is propagated using the state equation: $\tilde{\mathbf{x}}_{k}[i]=\mathbf{h}(\tilde{\mathbf{x}}_{k}[i-1],\underline{\mathbf{z}}[i];\boldsymbol{\theta})$ to simulate the evolution of this auxiliary state variable for $i=k-K,\cdots,k$. The prediction is then computed using $\hat{\mathbf{y}}[k]=\mathbf{g}(\tilde{\mathbf{x}}_{k}[k],\underline{\mathbf{z}}[k];\boldsymbol{\theta})$. The parameters are obtained by minimizing a cost function similar to (1).

Multiple shooting is equivalent, in the sense of Theorem 2 and Corollary 3, to solving the original (single shooting) problem regardless of the choice of simulation interval $\Delta m_{\text{max}}$. MSA-PEM, on the other hand, is equivalent to the original formulation only if $K=N$. The next example illustrates how the statistical properties of the method change as $K$ varies from $1$ to $N$.

6.2 Example: estimating the output error model for second-order under-damped system

A dataset with $N=300$ samples is generated using the equation:

[TABLE]

Here $\bar{y}$ represents the noiseless output, and $v$ represents the white output noise that is introduced during the data generation. At each $k$, $v[k]$ is a Gaussian random variable with standard deviation $\sigma_{v}=0.05$.

Let $\boldsymbol{\theta}=(\theta_{1},\theta_{2},\theta_{3})$ denote the parameter vector used for the data generation. Three different settings are considered: (a) $\boldsymbol{\theta}=(0.5,-0.2,2)$; (b) $\boldsymbol{\theta}=(1.5,-0.7,0.5)$; and, (c) $\boldsymbol{\theta}=(1.8,-0.95,0.1)$. The three settings correspond to underdamped linear systems with different response-times: in (a) the system responds fast to input changes and has a time constant $\tau=0.25$; in (c), it responds slowly and $\tau=0.9$; and, (b) is an intermediate setting with $\tau=0.75$.

From the synthetically generated data, the parameters of a model with the same structure are estimated. The experiment is repeated 100 times, each time corresponding to a different realization (i.e., different random seed) of the data generation. The result of the Monte Carlo procedure is displayed in Fig. 5 for the parameter estimation using ARX, standard (single shooting) output error, multiple shooting, and MSA-PEM. We show the difference between estimated and true values for the first parameter, $\hat{\theta_{1}}-\theta_{1}$. Similar results have been obtained for $\hat{\theta}_{2}$ and $\hat{\theta}_{3}$.

Due to the presence of white output error, using output error estimation yields consistent and well-behaved results. ARX estimation, on the other hand, is biased, and its histogram is not centered around zero. This is more evident either in 5(b) or in 5(c). In Fig. 5(a), which correspond to the fast response system, there is not much difference between the two estimation procedures, mostly because any disturbance fade away very quickly, yielding similar statistical properties for the two estimators.

Multiple shooting estimation results are similar to the (single shooting) output error estimation. For some choices of $\Delta\rm{m}_{\max}$, the optimization problem might actually be made harder due to the increase in the problem dimension. In this case, the algorithm seems to fail a few times, which can be observed in the histogram as a few outliers in the estimated parameter distribution. Except for these outliers, all the different choices of $\Delta\rm{m}_{\max}$ yield very similar distributions for the estimated parameters.

For MSA-PEM, the estimation results vary qualitatively among the three settings. For the intermediate setting, midway qualitative behavior between ARX and output error estimation is obtained, approaching that of output error estimation as $K$ increases. For the fast-response system, a similar interpretation is possible, however, there is not much difference between using ARX and output error estimation to begin with (cf. Fig 5 (a)), hence all $K$-step-ahead choices yield similar results. Finally, for the slow-response setting, minimizing the multi-step-ahead prediction error does not offer a direct compromise between ARX and output error estimation and, for $K=10$ or $K=20$, it is possible to observe a bimodal distribution with one peak distant from zero. And $K=7$ yields the best results.

In Fig. 6, the computational cost of multiple shooting and MSA-PEM are compared. The average running time divided by the number of function evaluations is displayed in Fig. 6(a) and (b). This accounts for the computational cost of computing the cost function, derivatives, and performing the matrix factorizations needed by the optimization algorithm at each iteration. Figure 6(c) and (d) display the number of cost function evaluations needed for the optimization algorithm to converge. The next section studies time and memory complexity of the methods and explains the obtained results in a more general context.

6.3 Computational cost

As observed in Fig. 6(b), MSA-PEM has an average running time per iteration that grows linearly with the number of stages. Since the problem dimension and the factorization cost per iteration remain constant, it is the time complexity of computing the cost function $V$ (and its derivatives) that accounts for the linear growth. In Appendix E, we discuss in detail the computation of $V$ and arrive at the time complexity of $\mathcal{O}(NK(n_{y}+n_{u}))$.

Multiple shooting estimation, on the other hand, has an average running time per iteration that decreases with the maximum simulation length $\Delta{\rm m}_{\max}$. This is displayed in Fig. 6(a). In this case, the cost of computing $V$ is $\mathcal{O}(N(n_{y}+n_{u}))$ and does not depend on the simulation length. The increase in the problem dimension, however, yields more costly factorizations (cf. Appendix A), which results in the increased computational cost per iteration.

For multiple shooting, shorter simulation lengths $\Delta{\rm m}_{\max}$ allow the solver to converge with fewer function evaluations. This is in agreement with the theory presented in this paper, which suggests longer simulation lengths might result in poor smoothness properties and make the optimization problem harder to solve. Applying the same reasoning to MSA-PEM, it is natural to expect that smaller values of $K$ would yield convergence with fewer iterations. This, however, does not seem to be the case for the slow time-response setting, maybe because MSA-PEM often converges to the wrong parameters in this setting (as shown in Fig. 5(i)).

In general, the number of function evaluations for MSA-PEM is lower than for multiple shooting because the resulting problem is unconstrained, rather than constrained, for which more efficient procedures are available.

Finally, MSA-PEM and multiple shooting are both amenable to parallelization, even though this is not explored in the examples. During MSA-PEM implementation, it is possible to instantiate $N$ different processes (or threads) for computing each prediction $\hat{y}[k]$ (and the corresponding derivatives) in parallel. Multiple shooting, on the other hand, subdivides the problem into $(N/\Delta{\rm m}_{\max})$ independent subproblems.

6.4 Cost function smoothness and estimation tradeoffs

Both MSA-PEM and multiple shooting allow the user to control the maximum simulation length, hence the smoothness of the cost function (cf. Section 3), by setting $K$ for the MSA-PEM and $\Delta{\rm m}_{\max}$ for the multiple shooting. These parameters, however, affect diverse aspects of the estimation problem, which are summarized in Table 2. Multiple shooting has an exact equivalence with the original single shooting problem regardless of the simulation length and the parameter $\Delta{\rm m}_{\max}$ offers a tradeoff between the problem dimension and smoothness properties of the cost function. MSA-PEM keeps the problem dimension fixed, and, as $K$ increases, it is possible to approach the behavior of methods with distinct statistical properties (not always in a smooth way, as shown in Fig. 5(i)), at the expense of an increased computational cost and worse smoothness properties.

Some refinements are proposed in the literature. For instance, [19] propose to average over different values of $K$, to obtain a slower but smoother convergence to output error estimation as $K\rightarrow N$. In [21], on the other hand, the method starts from the one-step-ahead prediction solution and increases $K$ by one in each iteration until convergence. We give an numerical example using such approach in Appendix F.

7 Conclusion

The relevance of this paper lies in the rather general setting for which the proposed methods and results hold. The main technical contribution is to show that for dynamic prediction models that are non-contractive (i.e. do not converge asymptotically to a single stable point) in the region of interest, the upper bound for the Lipschitz constant and the $\beta$-smoothness blows up exponentially with the simulation length, and this can make the optimization problem very hard to solve. This was illustrated using numerical examples with systems that are non-contractive due to the presence of chaotic regions and unstable equilibrium points. Because of these regimes, the objective function becomes very intricate in some regions of the parameter space and the optimization algorithm fails to find a good solution. Even for problems that are contractive in the region of interest, multiple shooting might help in preventing the solver from getting stuck in undesirable regions of the parameter space, hence facilitating the convergence to a good solution (cf. Section 5.2).

Multiple shooting makes the simulation length a design parameter and hence allows solving optimization problems that would be infeasible in a single shooting setting. The price paid compared to single shooting methods is that a nonlinear constrained optimization problem must be solved instead of an unconstrained one. It also makes it harder to generalize to situations other than batch training, such as online training. MSA-PEM is another approach that allows some control over the smoothness of the cost function, but with different tradeoffs (see Table 2) between smoothness, asymptotical properties and computational cost to be taken into consideration.

Appendix A Computing the derivatives

A.1 Sensitivity equations

Let the Jacobian matrices of $\mathbf{h}(\mathbf{x},\underline{\mathbf{z}};\boldsymbol{\theta})$ with respect to $\mathbf{x}$ and to $\boldsymbol{\theta}$ evaluated at the point $(\mathbf{x}[k],\underline{\mathbf{z}}[k];\boldsymbol{\theta})$ be denoted, respectively, as $A_{k}$ and $B_{k}$. Similarly, the Jacobian matrices of $\mathbf{g}(\mathbf{x},\underline{\mathbf{z}};\boldsymbol{\theta})$ are denoted as $C_{k}$ and $F_{k}$. Also, we denote the Jacobian matrices of $\hat{\mathbf{y}}[k]$ with respect to $\boldsymbol{\theta}$ and to $\mathbf{x}_{0}$ as $J_{\boldsymbol{\theta}}[k]$ and $J_{\mathbf{x}_{0}}[k]$. And the Jacobian matrices of $\mathbf{x}[k]$ are denoted as $D_{\boldsymbol{\theta}}[k]$ and $D_{\mathbf{x}_{0}}[k]$.

A direct application of the chain rule to (6) gives a recursive formula for computing the derivatives of the predicted output in relation to the parameters in the interval ${1\leq k\leq N}$:

[TABLE]

A similar recursive formula may be used for computing the derivatives of the predicted output in relation to the initial conditions:

[TABLE]

Finally, we define $D[k]=[D_{\boldsymbol{\theta}}[k],D_{\mathbf{x}_{0}}[k]]$ and

$J[k]=[J_{\boldsymbol{\theta}}[k],J_{\mathbf{x}_{0}}[k]]$.

A.2 Single shooting

For the cost function $V$ defined as in (1), its gradient $\nabla V$ is given by:

[TABLE]

Its Hessian $\nabla^{2}V$ is given by:

[TABLE]

where $\mathbf{S}[k]=\sum_{j=1}^{N_{y}}\hat{y}_{j}[k]\nabla^{2}\hat{y}_{j}[k]$. Ignoring $\mathbf{S}[k]$ is a common approximation used in least-squares algorithms, that will also be used here when computing derivatives numerically.

A.3 Multiple shooting

In order to solve the problem using the sequential quadratic programming solver [28] we must be able to compute: i) The cost function $V^{M}$; ii) its gradient $\nabla V^{M}$; iii) the constraints; iv) the Jacobian matrix of the constraints (which can be represented using a sparse representation); and, v) for any given vector $\mathbf{p}$, the product of the Lagrangian444The Lagrangian is given by: $\mathcal{L}(\boldsymbol{\phi},\boldsymbol{\lambda})=V(\boldsymbol{\phi})+\boldsymbol{\lambda}^{T}\mathbf{c}(\boldsymbol{\phi})$. Hessian and the vector $\mathbf{p}$ (the full Lagrangian Hessian matrix does not need to be computed). The following sequence provides a way of computing all derivatives required by the optimizer.

Algorithm 1** (Derivatives).**

For a given parameter $\boldsymbol{\theta}$ and set of initial conditions:

For $i=1,\cdots,M$, do:

1. (a)

For $k=\text{m}_{i}+1,\cdots,\text{m}_{i+1}$:

1. i.

Compute $\mathbf{x}^{i}[k]$ and $\hat{\mathbf{y}}^{i}[k]$ with (6). 2. ii.

Compute $A_{k}$, $B_{k}$, $C_{k}$ and $F_{k}$. 3. iii.

Compute $D^{i}[k]$ and $J^{i}[k]$ with the sensitivity equations. 2. (b)

Compute $V_{i}$ using Eq. (11). 3. (c)

Compute $\nabla V_{i}$ using a formula equivalent to (23). 4. (d)

Approximate the product of the Hessian with a given vector, $\nabla^{2}V_{i}\mathbf{p}$, using the first terms from a expression equivalent to (24) 2. 2.

Compute $V^{M}$ with (12); 3. 3.

Compute $\nabla V^{M}=\sum_{i=1}^{M}\tfrac{\Delta{\rm m}_{i}}{N}~{}\nabla V_{i}$; 4. 4.

Compute the value of the constraint from the values of $\mathbf{x}^{i}[\text{m}_{i+1}]$, $i=1,\cdots,M$; 5. 5.

Compute the Jacobian matrix of the constraints from $J^{i}[\text{m}_{i+1}]$, $i=1,\cdots,M$; 6. 6.

Compute $\nabla^{2}V^{M}\mathbf{p}=\sum_{i=1}^{M}\tfrac{\Delta{\rm m}_{i}}{N}\nabla^{2}V_{i}\mathbf{p}$; 7. 7.

Compute the product of the Hessian $\boldsymbol{\lambda}^{T}\mathbf{c}(\boldsymbol{\phi})$ with a vector $\mathbf{p}$ using 2-point finite differences; 8. 8.

Compute the product of the Lagrangian Hessian and a vector $\nabla^{2}\mathcal{L}(\boldsymbol{\phi},\boldsymbol{\lambda})\mathbf{p}$, summing the Hessians computed in steps 6 and 7.

Some approximations were used for computing the second derivatives:

1. for computing the Hessian of the objective function, the standard least-squares approximation for the Hessian is used; and, 2) for computing the Hessian of the constraint we use finite-difference approximation. The use of finite differences here comes inexpensively because we only need to evaluate the Hessian times a vector, and not the full matrix. Hence, it can be done at the cost of an extra Jacobian matrix evaluation.

Notice that step (1) from the above algorithm can be parallelized, with different processes (or threads) performing the computation for different values of $i$.

Appendix B Proofs

B.1 Preliminary results

Lemma 4**.**

Let $\mathbf{f}$ and $\mathbf{g}$ be two Lipschitz functions on $\Omega$ with Lipschitz constants $L_{f}$ and $L_{g}$. Then,

1. a)

$\mathbf{f}+\mathbf{g}$* is also a Lipschitz function on $\Omega$ with the (best) Lipschitz constant upper bounded by $(L_{f}+L_{g})$;* 2. b)

if, additionally, $\mathbf{f}$ and $\mathbf{g}$ are bounded by $M_{f}$ and $M_{g}$ on $\Omega$, then $\mathbf{f}\mathbf{g}$ is also a Lipschitz function on $\Omega$ with the (best) Lipschitz constant upper bounded by $(L_{f}M_{g}+L_{g}M_{f})$.

B.2 Proof of Theorem 1 (a)

Assume two different trajectories resulting from simulating the system (6) with parameters and initial conditions $(\mathbf{x}_{0},\boldsymbol{\theta})$ and $(\mathbf{w}_{0},\boldsymbol{\phi})$, respectively. We denote the corresponding trajectories by $\mathbf{x}$ and $\mathbf{w}$. Let us call:

[TABLE]

Because h and g are Lipschitz in ${(\mathbf{x},\boldsymbol{\theta})}$ we have:

[TABLE]

for all $(\mathbf{x},\underline{\mathbf{z}},\boldsymbol{\theta})$ and $(\mathbf{w},\underline{\mathbf{z}},\boldsymbol{\phi})$ in $(\Omega_{\mathbf{x}},\Omega_{\underline{\mathbf{z}}},\Omega_{\boldsymbol{\theta}})$. Applying these relations recursively we get that:

[TABLE]

Since $L_{h}$ is positive, the constant multiplying the second term in the above equation is always larger than the constant multiplying the first one. Hence, taking the square root on both sides of the above inequality and after simple manipulations, we get:

[TABLE]

where:

[TABLE]

Since $\Omega$ is compact and $\hat{\mathbf{y}}[k]$ is a (Lipschitz) continuous function of the parameters and initial conditions, then $\hat{\mathbf{y}}[k]$ is bounded in $\Omega$, i.e. $\|\hat{\mathbf{y}}[k]\|\leq M(k)$. And, it follows from (26) and from the existence of an invariant set555There are multiple ways to guarantee the invariant set premise will hold, but a very simple way is to just choose h such that $\textbf{h}(\mathbf{0},\underline{\mathbf{z}};\mathbf{0})=\textbf{0}$. In this case, $\{\mathbf{0}\}$ is an invariant set and if $\Omega_{\boldsymbol{\theta}}$ contain this point the premise is satisfied. For this specific case, one can just choose $[\boldsymbol{\phi},\boldsymbol{w}_{0}]=\mathbf{0}$ and it follows from (26) that $\|\hat{\mathbf{y}}[k]\|\leq L_{g}S(k)\|[\boldsymbol{\theta},\boldsymbol{x}_{0}]\|=\mathcal{O}(S(k))$. The more general case, for any invariant set, follows from a similar deduction. in $\Omega$ that $M(k)=\mathcal{O}(S(k))$.

The following inequality follows from (1):

[TABLE]

where ${L_{y}=\max_{1\leq k\leq N}\|\mathbf{y}[k]\|}$. And, by putting together (28) and (26):

[TABLE]

for ${L_{V}=\left(\tfrac{2L_{g}}{N}\sum_{k=1}^{N}(L_{y}+M(k))S(k)\right)}$. The asymptotic analysis of this expression with regard to $N$ yields (8).

B.3 Proof of Theorem 1 (b)

It follows from (23) that:

[TABLE]

where we have used the notation $\Delta J[k]$ to denote the difference between $J[k]$ evaluated at $(\boldsymbol{\theta},\boldsymbol{x}_{0})$ and $(\boldsymbol{\phi},\boldsymbol{w}_{0})$. Analogously, $\Delta(J[k]\hat{\mathbf{y}}[k])$ denote the difference between $J[k]\hat{\mathbf{y}}[k]$ evaluated at the two distinct points.

From equation (21) and (22) it follows that:

[TABLE]

Since, the Jacobian of h is Lipschitz with Lipschitz constant $L_{h}^{\prime}$, it follows that:

[TABLE]

Using a procedure analogous to the one used to get Eq. (26), it follows that:

[TABLE]

where $S(j)$ is defined as in (27). An identical formula holds for $B_{j}$ and a similar formula, replacing $L_{h}^{\prime}$ with $L_{g}^{\prime}$, holds for $C_{j}$ and $F_{j}$.

Since h and g are Lipschitz with Lipschitz constants $L_{h}$ and $L_{g}$ it follows that $\|A_{j}\|\leq L_{h}$, $\|B_{j}\|\leq L_{h}$, $\|C_{j}\|\leq L_{g}$ and $\|F_{j}\|\leq L_{g}$. Hence, it follows from (26), (30), (32) and the repetitive application of Lemma 4 that $\|\Delta J_{\boldsymbol{\theta}}[k]\|$, $\|\Delta J_{\mathbf{x}_{0}}[k]\|$, $\|\Delta(J_{\boldsymbol{\theta}}[k]\hat{\mathbf{y}}[k])\|$ and $\|\Delta(J_{\mathbf{x}_{0}}[k]\hat{\mathbf{y}}[k])\|$ are upper bounded by $\|[\boldsymbol{\theta},\boldsymbol{x}_{0}]^{T}-[\boldsymbol{\phi},\boldsymbol{w}_{0}]^{T}\|$ multiplied by the following constants:

[TABLE]

where $T(k)=(L_{g}^{\prime}M(k)+L_{g}^{2})$ and:

[TABLE]

Hence,

[TABLE]

where

[TABLE]

Putting everything together the asymptotic analysis of $L_{V}^{\prime}$ results in (9).

Appendix C Lipschitz analysis for the multiple shooting

The next theorem follows from basic inequality manipulation and relates the Lipschitzness and $\beta$-smoothness of the cost function $V^{M}$ with that of its components $V_{i}$.

Theorem 5**.**

Defining $V^{M}$ as in (12), if each component $V_{i}$ is Lipschitz continuous with constant $L_{V_{i}}$ then $V^{M}$ is also Lipschitz with constant equal to or smaller than $L_{V^{M}}=\max(L_{V_{1}},\cdots,L_{V_{M}})$. Additionally, if the gradient of each component $\nabla V_{i}$ is Lipschitz continuous with constant $L_{V_{i}}^{\prime}$ then $\nabla V^{M}$ is also Lipschitz with constant equal to or smaller than $L_{V^{M}}^{\prime}=\max(L_{V_{1}}^{\prime},\cdots,L_{V_{M}}^{\prime})$.

{pf}

For $\boldsymbol{\theta}_{\text{ext}}=(\boldsymbol{\theta},\mathbf{x}_{0}^{1},\cdots,\mathbf{x}_{0}^{M})$ and $\boldsymbol{\phi}_{\text{ext}}=(\boldsymbol{\phi},\mathbf{w}_{0}^{1},\cdots,\mathbf{w}_{0}^{M})$ we have that:

[TABLE]

where $L_{V^{M}}=\max(L_{V_{1}},\cdots,L_{V_{M}})$. And similarly, ${L_{V^{M}}^{\prime}=\max(L_{V_{1}}^{\prime},\cdots,L_{V_{M}}^{\prime})}$, which yields the second result.

Appendix D Asymptotic properties of prediction error methods

D.1 Notation and data generation process

Consider $\mathbf{y}[k]$ and $\mathbf{u}[k]$ to be one realization of the random variables $\mathbf{Y}[k]$ and $\mathbf{U}[k]$. And denote:

[TABLE]

Hence, rewriting the data generation difference equations for nonlinear ARX, output error and ARMAX (see Sections 2.1, 2.2 and 2.3) with the random variables yield, respectively:

[TABLE]

Here $\mathbf{V}[k]$ is a random variable representing the noise that is injected in the system. Notice that there is a deterministic additive relation between $\mathbf{Y}[k]$ and $\mathbf{V}[k]$. Hence, if $\mathbf{Y}[k]$ and $\mathbf{U}[k]$ are determined, so is $\mathbf{V}[k]$, or, conversely, if $\mathbf{V}[k]$ and $\mathbf{U}[k]$ are determined, so is $\mathbf{Y}[k]$.

D.2 Optimal output prediction

Let us define the optimal output prediction at time $k$ as the following conditional expectation:

[TABLE]

which is, in the least square sense, the best prediction for the output given its previous values.666 This prediction provides the smallest squared conditional expected error between the predicted and observed values: \hat{\mathbf{y}}_{*}[k]=\arg_{\hat{\mathbf{Y}}}\min E\left\{\|\mathbf{Y}[k]-\hat{\mathbf{Y}}\|^{2}~{}\Big{|}~{}\underline{\mathbf{u}}[k],\underline{\mathbf{y}}[k-1]\right\}.

For the nonlinear ARX, output error and ARMAX models the optimal output predictions are given, respectively, by:

[TABLE]

which follows from a direct application of the definition (33) to the stochastic difference equations that are assumed for the data generation in each case.

D.3 Ideal cost function

Ideally, the model predicted output $\hat{\mathbf{y}}[k]$ should be as close as possible to the optimal one $\hat{\mathbf{y}}_{*}[k]$. The distance between the two series is quantified by the cost function:

[TABLE]

Now, back to the examples, given a parametrized function $\mathbf{f}_{\boldsymbol{\theta}}$ with the correct model structure (that is: there exist $\boldsymbol{\theta}^{*}$ such that $\mathbf{f}_{\boldsymbol{\theta}^{*}}=\mathbf{f}^{*}$), it follows from Eqs. (3), (4) and (5) that $\boldsymbol{\theta}^{*}$ yields $V_{*}=0$, and hence $\boldsymbol{\theta}^{*}$ is a minimizer of $V_{*}$.

D.4 Uniform convergence of $V$

The optimal cost function $V_{*}$ is not available for optimization. Nevertheless, under some mild regularity conditions, it has been proved that $V\rightarrow V_{*}$ in probability as $N\rightarrow\infty$ and that this convergence is uniform [2].

If $V_{*}$ has a single minimum $\boldsymbol{\theta}^{*}$ and $V$ is convex in a convex set containing $\boldsymbol{\theta}^{*}$ the minimizer of $V$ converges to the minimizer of $V_{*}$ [31, Theorem 2.1]. Alternative conditions for this to hold are given in [2]. For our three examples, this would imply $\mathbf{f}_{\theta}\rightarrow\mathbf{f}^{*}$ (convergence in probability for $N\rightarrow\infty$). If the conditions are not satisfied, the uniform convergence, at least, guarantees that the minimizer has an equivalent performance (as $N\rightarrow\infty$). Additionally, in [32] it is shown that if the solution is unique the estimator has an asymptotic normal distribution.

Appendix E Computational cost of MSA-PEM

As mentioned in Section 6.3, the computation of the cost function $V$ (and its derivatives) for the MSA-PEM method has linear time complexity with $K$. The rate at which the time complexity grows, however, depends on implementation choices. We use the model estimation described in the experiment in Section 6.2 to explain the tradeoffs of those choices. Propagating $K=3$ steps-ahead the linear system $y[k]=\theta_{1}y[k]+\theta_{2}y[k-2]+\theta_{3}u[k-1]$ can be performed in two ways. The first is to introduce intermediate variables:

[TABLE]

The other is to simplify the computation in advance, resulting in a system of higher order:

[TABLE]

where the $\boldsymbol{\theta}^{\prime}$ coefficients depend on $\boldsymbol{\theta}$ (for instance $\theta_{1}^{\prime}=\theta_{1}^{3}-2\theta_{1}\theta_{2}$) and are computed in advance. For the first approach, computing the cost function has time-complexity $\mathcal{O}(NK(n_{y}+n_{u}))$; for the second, the computational cost is slightly smaller: $\mathcal{O}(N(K+n_{u}+n_{y}))$. The reduced computational cost is obtained because some computations are performed in advance and stored. Hence, the additional time-efficiency comes at the cost of storing the computation in advance, and the memory of the first method is $\mathcal{O}(n_{y}+n_{u})$, while in the second approach the memory complexity is $\mathcal{O}(K+n_{y}+n_{u})$.

For linear systems, it is easy to find the optimal simplified computation in advance by using Diophantine equations [19]. In the experiment in Section 6.2, however, we use the less time-efficient approach. The general tone of our discussion justifies this choice since this implementation can easily be applied to generic systems. And, while heuristic solutions might be employed to simplify the computations in advance, there is no general approach for all types of nonlinear systems. Regardless of this design choice, however, the observed linear growth of the computational cost with $K$ is expected.

Appendix F Additional experiments

F.1 Neural network for modeling pilot plant

This example makes use of data from the level process station described in Example 1 from [3]. As in the original paper, we use a neural network to model the water column height as a function of the voltage applied to a control valve that modulates the water flow. We compare three different training methods: i) minimizing the one-step-ahead prediction error (NN ARX) ii) minimizing the free-run simulation error using single shooting method (NN OE - SS); and, iii) minimizing the free-run simulation error using multiple shooting method (NN OE - MS).

The neural network (NN) training depends on the weight initialization, hence the performance of the neural network can be regarded as a random variable and is displayed in Fig. 7, which compares the empirical cumulative distribution of the mean square error (MSE) over the validation dataset for the three methods. A linear ARX model ($n_{y}=1$ and $n_{u}=1$) was trained and tested under the same conditions to serve as the baseline. Methods (i) and (ii) and the linear ARX baseline were described in [3]. Method (iii) is introduced in this paper.

The cumulative distribution function gives, for each $x$-axis value, the probability of the method to yield a validation MSE smaller than or equal to this value. It was estimated from 100 realizations of the neural network training procedure. Fig. 7 shows that for more than 90% of the realizations, estimating the parameters by minimizing the free-run simulation (NN OE) offers significant advantages over the minimization of the one-step-ahead error (NN ARX). When using a standard single shooting formulation, however, it also makes the parameter estimation procedure more sensitive to the initial conditions, with the algorithm yielding some really bad results for some initial choices [3]. This results in a long-tailed distribution for the MSE (Fig. 7). More precisely, in 10 out of 100 realizations the NN OE - SS model yields a performance that is inferior to the linear ARX baseline, some of the realizations worse than the linear baseline by a factor of 100. The performance of the NN OE - SS and NN OE - MS is very similar for 90% of the realizations, the tail of the distribution, however, is very different, with the multiple shooting procedure rarely producing very bad results. In order to highlight the differences, results where NN OE - SS and NN OE - MS are worse than the baseline are presented, respectively, as blue and red circles in Fig. 7.

This example illustrates how the use of multiple shooting alleviates the problem of high sensitivity to initial conditions, making it possible to estimate output error models with extra robustness against variations of the initial conditions and lower probability of getting trapped at local minima with very bad performance.

This example also shows the limitations of the multiple shooting formulation. The training time for NN OE - MS model is 282 seconds, for NN OE - SS model is 3.9 seconds, and for NN ARX model is 3.3 seconds. This means that the single shooting parameter estimation could be repeated, roughly, 70 times for each multiple shooting run. Hence solving the single shooting problem several times and choosing the best result would also avoid very bad solutions and could, still, be computationally less expensive than solving the multiple shooting problem. The longer training time is due to two factors: i) per iteration the multiple shooting approach takes, roughly, 3.5 times more than the single shooting approach; and, ii) it takes, approximately, 20 times more iterations to converge. Both are consequences of the fact that a higher dimensional constrained optimization problem is being solved.

F.2 Incremental initialization approach for the MSA-PEM

Here, we revisit the setup of the first numerical example presented in [21]. A dataset with $N=500$ samples is generated using the equation:

[TABLE]

The data generator parameters are $\theta_{1}=0.6$ and $\theta_{2}=-0.5$, $\bar{y}$ represents the noiseless output and $v$ the white output noise. The input $u$ is generated by the following AR process:

[TABLE]

where $\eta$ is white Gaussian noise with variance one. The output noise $v$ has variance $0.25\lambda_{y}^{2}$, where $\lambda_{y}^{2}=0.18^{2}$ is the output variance.

A model with the same structure as the data generator is estimated using MSA-PEM. Fig. 8 displays the cost function and the result of parameter estimation for several initial optimization conditions. Besides the results of the vanilla implementation, in blue, we also show, in black, the refinement proposed in [21]. For this refinement, the prediction starts from the one-step-ahead solution, and in each iteration increases $K$ by one. That is, the $K$ MSA-PEM problem is initialized from the (partial) solution of the $K-1$ MSA-PEM estimation. The intermediary results for the refined implementation are displayed in the plot. With the algorithm being iterated up to the $K$-th iteration.

The effect of increasing the simulation length can be observed in this example, with a larger $K$ yielding less smooth cost functions. The refinement proposed in [21] allows the estimation to approach the desired asymptotic statistical properties while optimizing more complicated cost functions progressively. The results in Fig 8, show this refinement yield improvements over the vanilla implementation, maybe because it avoids optimizing the harder optimization problems from a distant initial estimate.

Appendix G Equality-constrained optimization solver

This appendix describes the implementation of the solver described in [28], which is used in the numerical examples. The method is able to solve large-scale equality-constrained optimization problems:

[TABLE]

where $V$: $\mathbb{R}^{n}\rightarrow\mathbb{R}$ and $\mathbf{c}$: $\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ are twice continuously differentiable functions.

The algorithm solves a sequence of Taylor approximations in order to gradually converge towards a local solution of (37). At the $k$-th iteration, the algorithm builds a local model around the current iterate $\boldsymbol{\theta}_{k}$, computes the update $\mathbf{p}_{k}$ by solving a quadratic programming problem and then updates the solution:

[TABLE]

The algorithm is a trust-region method [33], in the sense that the update must respect $\|\mathbf{p}_{k}\|\leq\Delta_{k}$, for which $\Delta_{k}$ is known as trust radius and should reflect the trust the algorithm has on the local approximation of the cost function. If the current local approximation is not a good one, the method does not allow a very large step to be taken.

G.1 Quadratic programming subproblem

At each iteration, in order to compute the step update $\mathbf{p}_{k}$, the method solves the trust-region quadratic programming (QP) sub-problem:

[TABLE]

where $\boldsymbol{\theta}_{k}$ is the current iterate, $\boldsymbol{\lambda}_{k}$ is the current approximation of Lagrange multipliers, $\nabla V(\boldsymbol{\theta}_{k})^{T}$ is the gradient of $V$ and $J(\boldsymbol{\theta}_{k})$ is the Jacobian matrix of $\mathbf{c}$. As in Appendix A, $\mathcal{L}(\boldsymbol{\theta}_{k},\boldsymbol{\lambda}_{k})^{T}$ denotes the Lagrangian:

[TABLE]

and $\nabla^{2}_{\theta\theta}\mathcal{L}(\boldsymbol{\theta}_{k},\boldsymbol{\lambda}_{k})^{T}$ is the Hessian in relation to the variable $\boldsymbol{\theta}$.

The current approximation of the Lagrange multipliers $\boldsymbol{\lambda}_{k}$ is obtained by solving the least-squares problem:

[TABLE]

The vector $\mathbf{r}_{k}$ is included in the problem in order to guarantee linear and trust-region constraints are always compatible. $\mathbf{r}_{k}$ is defined as $J(\boldsymbol{\theta}_{k})\mathbf{v}_{k}+\mathbf{c}(\boldsymbol{\theta}_{k})$ for a value of $\mathbf{v}_{k}$ satisfying:

[TABLE]

for $0<\eta<1$ (in our implementation $\eta=0.8$). For this choice of $\mathbf{r}_{k}$ the linear constraints are always compatible with trust-region constraints.

The sub-steps are solved in rather economical fashion using efficient methods to get an inexact solution to each of the sub-problems. The QP problem (39) is solved using a variation of the dogleg procedure (described in [34], p.886) and (42) is solved using the projected conjugate gradient (CG) algorithm [35].

G.2 Implementation details

Here we define the merit function, $\phi$, which combines the constraints and the objective function into a single number that can be used to compare two points and to reject or accept a given step. It is given by:

[TABLE]

where the penalty parameter $\mu$ is updated through the iterations. This parameter must increase monotonically for the algorithm to converge. Some additional guidelines are provided in [34], p.891.

The selection of the trust radius $\Delta_{k}$ and the step rejection mechanism are based on the ratio:

[TABLE]

between the reduction of the merit function and the reduction predicted by the local (Taylor approximation) model:

[TABLE]

The ratio measures the agreement between the observed and expected reduction and is used as decision variable when accepting or rejecting the update and enlarging or reducing the trust radius $\Delta_{k}$.

G.3 Algorithm Overview

The full algorithm is summarized next.

Algorithm 2** (Trust-region sequential QP solver).**

*At each iteration, until some stop criterion is met (e.g.

${\|\nabla_{x}\mathcal{L}(\boldsymbol{\theta}_{k},\boldsymbol{\lambda}_{k})\|_{\infty}<10^{-8}}$), repeat:*

Compute $V(\boldsymbol{\theta}_{k})$, $\nabla V(\boldsymbol{\theta}_{k})$, $c(\boldsymbol{\theta}_{k})$ and $J(\boldsymbol{\theta}_{k})$; 2. 2.

Compute least squares Lagrange multipliers $\boldsymbol{\lambda}_{k}$; 3. 3.

Compute $\nabla^{2}_{xx}\mathcal{L}(\boldsymbol{\theta}_{k},\boldsymbol{\lambda}_{k})$; 4. 4.

Apply dogleg method in order to compute $v_{k}$ and $r_{k}$ (such that the resulting problem is feasible); 5. 5.

Compute $p_{k}$ using the projected CG method; 6. 6.

Choose penalty parameter $\mu_{k}$; 7. 7.

Compute reduction ratio $\rho_{k}$; 8. 8.

Accept or reject step $p_{k}$ using $\rho_{k}$ as decision variable; 9. 9.

Enlarge or reduce trust-radius using $\rho_{k}$ as decision variable.

Code avaibility

The code for reproducing the examples is available in: github.com/antonior92/MultipleShootingPEM.jl.

Acknowledgements

This work has been supported by the Brazilian agencies CAPES - Coordenacão de Aperfeiçoamento de Pessoal de Nível Superior (Finance Code: 001), CNPq - Conselho Nacional de Desenvolvimento Científico e Tecnológico (contract number: 303412/2019-4, 200931/2018-0 and 142211/2018-4) and FAPEMIG - Fundação de Amparo à Pesquisa do Estado de Minas Gerais (contract number: TEC 1217/98), by the Swedish Research Council (VR) via the projects NewLEADS – New Directions in Learning Dynamical Systems (contract number: 621-2016-06079) and Learning flexible models for nonlinear dynamics (contract number: 2017-03807), and by the Swedish Foundation for Strategic Research (SSF) via the project ASSEMBLE (contract number: RIT15-0012).