1 Introduction↩︎

Mechanical, electrical, chemical and biological processes are often dynamic with highly nonlinear response to disturbance, actuation, and initial state. Regulation of such processes requires accurate models of this nonlinear behaviour. Ideally, the dynamics of the process can be modelled through the combination of well-established physical principles. In some cases, however (e.g. biological systems), the physical principles governing the behaviour are speculative, over-simplified, or only partially understood. Furthermore, even when the physics of the problem are well-understood, there may still be substantial uncertainty in parameters of the model due to natural variation, errors in measurement, oversimplification of the physics or difficulty in direct measurement of these parameters. In such cases, there may be only a few uncertain parameters in the model, and they may be constrained to lie in some set (e.g. positivity is a common form of parameter constraint).

When it is not possible to measure the parameters of the system directly, these parameters must be inferred indirectly by examination of the response of the system to variations in input and initial state [1]. In the extreme case, where there is no understanding of the physics, we have black-box system identification [2] — e.g. Koopman operators [3], [4]; extended dynamic mode decomposition (EDMD) [5]; sparse identification of non-linear dynamics (SINDy) [6]–[8]; autoregressive models [9]; and neural networks [10].

When the physics are only partially understood or the system parameters are constrained, however, we have grey-box system identification [11] – a more challenging problem than black-box modeling. Specifically, suppose the system to be identified is known to have the form \[\begin{align} \dot{x}(t)&= f(t, x(t), \theta) \qquad x(0) = x_0\notag \\ y(t) &= g(t, x(t), \theta) \label{eqn:nonlinearODE} \\[-8mm]\notag \end{align}\tag{1}\] where the functions \(f\) and \(g\) are given, \(x(t) \in \mathbb{R}^{n_x}\) is the internal state, \(x_0\) is the initial state, and \(y(t)\in \mathbb{R}^{n_y}\) is some measurable output. The vector, \(\theta\in \Theta \subset \mathbb{R}^{n_\theta}\), represents the unknown parameters in the model. The dependence of \(f\) and \(g\) on time, \(t\), is typically used to represent the effect of known inputs, \(u(t)\) — in which case we might equivalently write the functions as \(f(x(t), u(t), \theta)\) and \(g(x(t), u(t), \theta)\).

Of course, when measurements of \(\dot{x}(t)\) are available, or \(x(t)\) is measurable (\(g=f\) or \(g(t,x)=x\), respectively) and the sampling instances are sufficiently dense, it is relatively straightforward to use least squares or some other form of regression to estimate the parameters, \(\theta\) [1]. However, when sampling instances are sparse in time or only output measurements are available, the system identification problem becomes significantly more challenging [2].

Specifically, in this paper, we consider the case of grey-box parameter estimation based on noisy irregular sampling of system output, \(y(t)\), and where the data may be generated using multiple experiments, each with different initial condition \(x_{0,i}\) and sampling instances. Note that implicit in this framework is that the initial condition, \(x_{0,i}\) for each experiment is unknown and must be estimated along with the system parameters, \(\theta\) (similar to the inverse problem [12]).

For a given parameterized vector field, \(f(t,x,\theta)\), we suppose that Eqn. 1 is well-posed and define the corresponding solution map, \(\phi_{f}(t, x, \theta)\), as the unique function which satisfies \[\partial_t \phi_{f}(t, x, \theta) = f(t, \phi_{f}(t, x, \theta), \theta),\quad \phi_{f}(0, x, \theta) = x \] for all \(t \ge 0\), \(x \in \mathbb{R}^{n_x}\) and \(\theta\in \Theta\). While it is almost never possible to obtain an analytic expression for the solution map of a nonlinear system, the \(\phi_{f}\) notation allows us to efficiently represent the problem of data-based system identification. Specifically, suppose that experimental data is available from solutions of Eqn. 1 consisting of measurements \(y_i\) at a sequence of discrete times \(\{t_i\in\mathbb{R}\}_{i=1}^{N_s}\) so that \(y_i = (1+n_i) g(t_i, \phi_{f}(t_i, x, \theta), \theta) + m_i\) where \(\theta,x\) are unknown and \(m_i\in \mathbb{R}^{n_y}\) and \(n_i \in \mathbb{R}\) represent measurement errors. We denote the set of all such measurements as \(Y = \{(y_i, t_i) \in \mathbb{R}^{n_y}\times\mathbb{R}\; : \; i = 1, ..., N_s\}.\)

We may now formulate the parameter estimation/system identification problem as \[\min_{\substack{\theta\in\Theta, x\in\mathbb{R}^{n_x}}} L_{ls}(x, \theta, \phi_{f}(\cdot, x, \theta)), \label{eqn:optimization}\tag{2}\] where for simplicity we use the least-squares Loss Function [13] defined as \[L_{ls}(x, \theta, \phi_{f}(\cdot, x, \theta))=\frac{1}{2N_s} \sum_{i=1}^{N_s} \|y_i - g(t_i, \phi_{f}(t_i, x, \theta), \theta)\|_2^2. \]

Problem 2 is an unconstrained optimization problem [14] in parameters \(x,\theta\). However, for nonlinear systems, the inclusion of the nonlinear solution map, \(\phi_f\), typically makes the problem non-convex – implying that descent-based algorithms will obtain at best local optimality. Furthermore, because explicit expressions for the solution map \(\phi_f\) are rarely available, obtaining even locally optimal solutions to this problem is challenging.

One approach to solving Problem 2 is to use numerical simulation in place of the solution map [15] and then apply gradient-free, black box optimization tools such as numerical gradient-approximation [16], parameter cascade [17], [18], nonparametric estimators [19]–[21], genetic algorithms [22], [23], and simulated annealing [24]. Unfortunately, however, such methods are computationally inefficient and may fail to converge to even locally optimal solutions.

By contrast, gradient-based methods are guaranteed to converge to local optima but require some method for computing the gradient of the solution map, \(\nabla_{x, \theta} \phi_{f}\). Because analytic expressions for the solution map are rarely available, however, such a gradient descent approach is not practical.

Faced with the limitations of black-box optimization and lack of analytic solution maps, and motivated by the use of gradient descent to obtain local optima, we consider now an alternative formulation of the parameter estimation/system identification problem posed in Eqn. 2 , which does not require knowledge of the solution map. Specifically, in Sec. 2, we reformulate the unconstrained optimization problem as the equivalent constrained optimization problem defined as \[\begin{align} \label{eqn:intro95optimization95problem95constrained} \min_{\substack{x \in X, \theta \in \Theta, \mathbf{u} \in C(\mathbb{R})}}&\; L_{ls}(x, \theta, \mathbf{u})\\[-3mm] \text{s.t.} &\; {\mathbf{u}}(t) = x+\int_{0}^t f(s,\mathbf{u}(s),\theta)ds \qquad \forall t \ge 0. \notag\\[-8mm] \notag \end{align}\tag{3}\] In Optimization Problem 3 , we have replaced the solution map \(\phi_f\) by the variable \(\mathbf{u}\) and added the constraint that \(\mathbf{u}\) be a solution of the given system. Clearly, Problems 2 and 3 are equivalent. Moreover, the constraint, \({\mathbf{u}}(t) = x+\int_{0}^t f(s,\mathbf{u}(s),\theta)ds\), may be conveniently represented using the Picard operator as \(\mathbf{u}=\mathcal{P}_{x, \theta} \mathbf{u}\), where \[(\mathcal{P}_{ x, \theta} \mathbf{u})(t):= x + \int_{0}^t f(s, \mathbf{u}(s), \theta) ds.\] Although Optimization Problem 3 does not require the solution map, it introduces an infinite-dimensional variable \(\mathbf{u}\) and associated equality constraint which must hold for all times, \(t\), on which the data is defined. Solving such problems without an explicit representation of the feasible set is known to be hard [25].

To address these challenges, in Sec. 3 we propose a new class of gradient-contractive algorithms (Alg. [alg3:modified95picard]) for solving the formulation of the parameter estimation problem in Eqn. 3 . In this algorithm, each iteration has two steps – gradient computation and contraction to the feasible set. Specifically, the first step computes the gradient of the loss function and updates the initial condition and parameter variables as \(x_{k+1},\theta_{k+1}\) where \[\begin{align} x_{k+1} &= x_{k} - \alpha \nabla_{x} L_{ls}(x_k, \theta_k, \mathcal{P}_{ x_k, \theta_k}^n \mathbf{u}_{k})\\ \theta_{k+1} &= \theta_k - \alpha \nabla_{\theta} L_{ls}(x_k, \theta_k, \mathcal{P}_{ x_k, \theta_k}^n \mathbf{u}_{k}), \\[-8mm] \end{align}\] and where the \(n^{th}\) order Picard iteration, \(\mathcal{P}_{ x, \theta}^n\), is introduced into the objective so that the gradient accounts for multipliers corresponding to the KKT conditions for the constraint \(\mathbf{u}=\mathcal{P}_{x, \theta} \mathbf{u}\) and where the step size \(\alpha > 0\) is computed through a line search. The second step contracts the variable \(\mathbf{u}_{k}(t)\) to the new feasible set (determined by \(x_{k+1}\) and \(\theta_{k+1}\)) using the contractive property of the Picard iteration as \[\mathbf{u}_{k+1}(t) = (1- \sigma) \mathbf{u}_{k}(t) + \sigma (\mathcal{P}_{ x_{k+1}, \theta_{k+1}} \mathbf{u}_k)(t) \] for some \(\sigma \in (0, 1]\). This second step translates the solution closer to the infinite-dimensional manifold described by the feasible set of Problem 3 – imitating the projection step of gradient-projection algorithms and resolving the problems introduced by the infinite-dimensional nature of the equality constraints. This approach also eliminates the need to explicitly parameterize the infinite-dimensional variable, \(\mathbf{u}\).

A limitation of this approach, of course, is the need for the Picard operator to be contractive on the interval of time on which the data was sampled. However, as described in Sec. 3, an extended form of Picard operator may be used, which is contractive over arbitrary intervals (assuming well-posedness of the solution map), but requires the inclusion of additional decision variables for each such extension.

To establish convergence of the proposed class of algorithms to local optimality, we may impose regularity and convexity on the objective function and vector field. However, for nonlinear vector fields, the convexity condition will likely fail for \(n\ge 1\). Hence, we necessarily limit our analysis to local convergence to locally optimal solutions. Such a limitation is not particularly strong, however, since many parameter identification problems in nonlinear systems may not have a unique solution or such parameters may be locally but not globally identifiable [26].

Note that Picard iteration has previously been used for analysis of black-box system identification [27], and parameter estimation in [28], [29]. Specifically, in [28] gradient descent of the unconstrained Problem 2 was considered for signal transduction networks, wherein a Picard iterate was used in lieu of the solution map (See Alg. 1). Unfortunately, convergence of the Picard iteration is typically limited to short time intervals, thus limiting the generality of this approach. An attempt to extend this approach to longer time intervals using multiple shooting methods was proposed in [29]. However, convergence of neither method has been proven or extended to sparse data, irregular sampling times or measured outputs.

In contrast to the work in [28], [29], the use of a constrained formulation of the parameter estimation problem, the use of gradient contractive algorithms, and the proposed extended form of Picard iteration (Alg. [alg4:ext95picard]) used in this paper allow for establishment of convergence proofs, arbitrary sampling intervals, and measured outputs. To establish these results, in Sec. 4, we first consider the general class of two-step gradient contractive algorithms and establish conditions for convergence based on properties of the contractive operator, \(\mathcal{P}\). Next, in Sec. 5, we examine the Picard operator \(\mathcal{P}_{ x, \theta}\) and establish conditions for convergence to the solution map and its gradient. In Sec. 6, we show these results apply to the extended Picard operator (Alg. [alg4:ext95picard]) and provide conditions for convergence of Alg. [alg4:ext95picard] to a local solution of the parameter estimation problem. The results are then applied to a rigorous battery of numerical tests, including: Van der Pol oscillator to evaluate the effect of the regularization parameter; the FitzHugh-Nagumo Neuron to evaluate the effect of irregular sampling; the Rosenzweig-MacArthur predator-prey model to evaluate the effect of sparse data; a tumor growth model to investigate the question of identifiability; and the Lorenz system for comparison with SINDy and black-box system optimization methods.

Notation↩︎

We use \(\mathbb{N}\), \(\mathbb{R}\), and \(\mathbb{R}_{+}\), to denote the natural numbers, the real numbers, and the non-negative real numbers, respectively. For a given \(J \in \mathbb{N}\) we denote \(\overline{1, J} = \{j \in \mathbb{N}\;|\; 1 \leq j\leq J\}\). For \(a, b \in \mathbb{R}_+\) we denote \({\lfloor a/b \rfloor} :=\underset{m\in\mathbb{N}, a \geq mb}{\texttt{argmax}} m\) and \(a \bmod{b} := a- \lfloor a/b\rfloor b\). For a given sets \(X, Y\) we denote \(\mathcal{F}(X, Y)\) to be the set of all functions \(f:X\rightarrow Y\). For compact \(S \subset \mathbb{R}^n\), we denote \(C^k(S)\) to be the space of \(k\) times continuously differentiable functions, \(f \in \mathcal{F}(S,\mathbb{R}^m)\) with norm \(\|f\| = \sup_{s\in S} \|f(s)\|_2\). We use \(C_{a}(S)\subset C(S) := C^0(S)\) to denote the ball of radius \(a\) as \(C_a(S) = \{\mathbf{u}\in C(S) \; |\; \|\mathbf{u}\| \leq a\}\) and use \(C_a\) when the domain is clear from context. For \(u \in C^1(S)\), we use \(\nabla_x u\) to denote the gradient as \(\nabla_x u = \begin{bmatrix}\frac{\partial}{\partial x_1}u, \dots, \frac{\partial}{\partial x_m}u \end{bmatrix}^T\). For compact \(X \subseteq \mathbb{R}^{n}\) we denote \(\Pi_X:\mathbb{R}^{n} \rightarrow X\) to be the projection operator so that \(\Pi_X(y) = \underset{x \in X}{\texttt{argmin}} \|x- y\|_2\). We say \(f\in \mathcal{F}(X,\mathbb{R})\) is strongly convex on \(X\subset \mathbb{R}^n\) with modulus \(\mu > 0\) if \(f(x) - \frac{\mu}{2} \|x\|^2\) is convex on \(X\). We say \(f\in \mathcal{F}(X\times Y,\mathbb{R})\) is strongly convex on \(X\subset \mathbb{R}^n\) uniformly in \(y\) if there exists \(\mu\geq 0\) such that \(f(x, y) - \frac{\mu}{2}\|x\|^2\) is convex on \(X\) for all fixed \(y \in Y\subset \mathbb{R}^m\). For \(X\subset \mathbb{R}^n\) and \(Y\subset \mathbb{R}^m\), we say that \(f\in \mathcal{F}(X , Y)\) is Lipschitz continuous with constant \(K\) if for all \(x_1, x_2 \in X\) we have \(\|f(x_1) - f(x_2)\| \leq K \|x_1 - x_2\|\) and for \(Z \subset \mathbb{R}^p\), \(f\in \mathcal{F}(X\times Z,Y)\) is Lipschitz continuous with respect to \(x\) with constant \(K\) if for all \(x_1, x_2 \in X\) and \(z \in Z\) we have \(\|f(x_1, z) - f(x_2, z)\| \leq K \|x_1 - x_2\|\).

2 Formulating the Parameter Estimation Problem for Nonlinear Differential Equations↩︎

Consider a parameterized nonlinear differential equation model in state-space representation. \[\begin{align} \dot{\mathbf{u}}(t) &= f(t, {\mathbf{u}}(t), \theta) \qquad {\mathbf{u}}(0) = x \notag\\ y(t) & = g(t, {\mathbf{u}}(t), \theta) \label{eqn:ODE} \\[-7mm] \notag \end{align}\tag{4}\] where \(t\in \Gamma\) is time (on interval \(\Gamma\)), \({\mathbf{u}}(t) \in \mathbb{R}^{n_x}\) is system state, \(\theta \in \Theta\) are unknown parameters (restricted to compact convex \(\Theta \subset \mathbb{R}^{n_\theta}\)), \(x\) are unknown initial conditions (restricted to compact convex \(X \subset \mathbb{R}^{n_x}\)), and \(y(t) \in \mathbb{R}^{n_y}\) are measured outputs. The function \(f\in C^1(\Gamma \times\mathbb{R}^{n_x}\times\Theta)\) is the vector field and \(g\in C^1(\Gamma\times\mathbb{R}^{n_x} \times\Theta)\) is the state-to-output map. We assume existence and uniqueness, and continuous dependence of solutions on initial conditions and parameters, implying existence of a solution map, \(\phi:\Gamma\times X\times \Theta \rightarrow \mathbb{R}^{n_x}\), which satisfies \[\partial_t \phi(t, x, \theta) = f(t, \phi(t, x, \theta), \theta), \quad \phi( 0, x, \theta ) = x. \label{eqn:solutionmap} \tag{5}\]

Now suppose we are given a set of measurements of the dynamical system for some initial condition \(x\) and true parameter value, \(\theta\) as \[Y =\{ (y_i, t_i) \; | \; y_i = (1+n_i)g(t_i, \phi(t_i, x, \theta), \theta) + m_i,\;i \in \overline{1, N_s} \},\] where \(m_i\in \mathbb{R}^{n_y}\) and \(n_i \in \mathbb{R}\) represent measurement errors.

The parameter estimation problem is to find initial condition, \(x \in X\), and parameter values, \(\theta\in \Theta\), that minimize the least squares error in the predicted output relative to the set of measurements. This can be formulated as either an unconstrained or a constrained optimization problem depending on whether we know the solution map. If the solution map, \(\phi\), is known, the problem can be formulated as an unconstrained optimization problem of the form \[\min_{\substack{x \in X, \theta \in \Theta}} L_{ls}(x, \theta, \phi(\cdot, x, \theta)), \label{eqn:optimization95problem95unconstrained}\tag{6}\] where \[L_{ls}(x, \theta, \phi(\cdot, x, \theta)) = \frac{1}{2N_s}\sum_{i=1}^{N_s} \|y_i - g(t_i, \phi(t_i, x, \theta), \theta)\|^2.\] Since the solution map is unknown, however, we may reformulate the parameter estimation problem by instead introducing a variable, \(\mathbf{u}\), representing the solution of the ODE and add the constraint that this variable be a solution of the ODE. \[\begin{align} \label{eqn:optimization95problem95constrained} \min_{\substack{x \in X, \theta \in \Theta, \mathbf{u} \in C(\Gamma)}}&\; L_{ls}(x, \theta, \mathbf{u})\\[-3mm] \text{s.t.} &\; {\mathbf{u}}(t) = x+\int_{0}^t f(s,\mathbf{u}(s),\theta)ds, \notag \\[-8mm]\notag \end{align}\tag{7}\] where the integral constraint is a necessary and sufficient condition for \(\mathbf{u}(t)\) to satisfy the ODE (Eqn. 4 ) with initial condition \(\mathbf{u}(0)=x\) and parameter values \(\theta\). The advantage of this formulation is that it does not require knowledge of the solution map. However, the disadvantage is that we have introduced an infinite-dimensional optimization variable in the form of \(\mathbf{u}\).

2.1 Reformulation using the Picard Operator↩︎

The unconstrained optimization problem in Eqn. 6 is expressed in terms of the solution map (which is unknown) and the constrained optimization problem in Eqn. 7 is expressed in terms of the variable \(\mathbf{u}\) (which is an infinite-dimensional variable and where the constraint must hold at an infinite number of times). The problems with both these formulations may be ameliorated through use of what is known as the Picard operator, \(\mathcal{P}_{x, \theta}\).

The critical property of the Picard operator is that if \(T\) is sufficiently small, iterations of the Picard operator converge to the solution of \(\dot{x}=f(x)\). That is, for any \(\mathbf{u}(t)\), \(\mathcal{P}_{x, \theta} \circ \cdots \circ \mathcal{P}_{x, \theta} \mathbf{u} \rightarrow \mathbf{u}^*\) where \(\mathcal{P}_{x, \theta} \mathbf{u}^*=\mathbf{u}^*\).

The Picard operator allows us to approximate Optimization Problem 6 for fixed \(n \in \mathbb{N}\) as \[\label{eqn:optimization95problem95unconstrained95approximated} \min_{\substack{x \in X, \theta \in \Theta}} L_{ls}(x, \theta, \mathcal{P}_{ x, \theta}^n \mathbf{u}), \tag{9}\] and to restate Optimization Problem 7 as \[\label{eqn:optimization95problem95constrained95approximated} \min_{\substack{x \in X, \theta \in \Theta, \\ \mathbf{u} \in C(\Gamma)}}L_{ls}(x, \theta, \mathbf{u})\quad \text{s.t.} \quad \mathcal{P}_{ x, \theta} \mathbf{u} = \mathbf{u}.\tag{10}\] Now, the approximation of unconstrained Optimization Problem 6 in Problem 9 has eliminated the need for knowledge of the solution map, but introduced a very complicated dependence on initial state, \(x\) and parameter, \(\theta\). Meanwhile, the restatement of constrained Optimization Problem 7 in 10 does not seem to have any immediate benefits. As will be discussed in the following section, however, the contractive properties of the solution map will allow us to eliminate the need for infinite-dimensional variables and constraints in 10 .

3 Gradient Descent for Parameter Estimation↩︎

In this section we consider algorithms to solve the optimization problems in Sec. 2. For the unconstrained problem (Problems 6 and 9 ), we focus on calculating the gradient of high-order Picard iterations. For the constrained optimization problem (Problems 7 and 10 ), we focus on an iteration which eliminates the infinite-dimensional nature of the constraint and variables.

3.1 A Gradient Descent Algorithm for the Unconstrained Formulation↩︎

As indicated in Subsec. 2.1, Optimization Problem 6 can be approximated using the Picard iteration as in 9 , where the Picard operator in this case is understood to act on the set of functions \(\mathbf{u}(t,x,\theta)\) and convergence of Picard iterates is then to the solution map. Unconstrained optimization problems of this form can be solved using gradient descent (combined with projection onto the convex set of allowable parameters, \(\theta\in \Theta\) and initial conditions, \(x \in X\)). Specifically, let \(\alpha > 0\) be a step size and recall that the projection operator \({\Pi}_X : \mathbb{R}^{n_x} \rightarrow \mathbb{R}^{n_x}\) is defined as \[\label{eqn:projection95operator} {\Pi}_X(z) := \underset{x \in X}{\texttt{argmin}} \|z - x\|^2_2\tag{11}\]

We may now propose Alg. 1, where the gradient is computed using the n-th order Picard iterate approximation of the solution map. As stated in the following theorem, this gradient also converges to the gradient of the solution map – implying that Alg. 1 is actually an approximation of the gradient descent algorithm as applied the original formulation of the problem in Eqn. 6 .

Convergence of the gradient descent in Alg. 1 to the solution of Optimization Problem 6 as \(n \rightarrow \infty\) is then guaranteed if sets \(X, \Theta\) are convex and the loss function is strongly convex – See, e.g., Thm. 4.32 in [30]. Of course, in reality, for most parameter estimation problems the dependence of the loss function in Eqn. 9 on the system parameters and initial state is not convex, and hence Alg. 1 is only guaranteed to converge to a set of local minima.

More significantly, convergence is only guaranteed in the limit \(n \rightarrow \infty\). However, when \(n\) becomes large, computation of the gradient \(\nabla_{x,\theta} {\mathcal{P}}_{x, \theta}^n\) becomes difficult. For this reason, we now consider an algorithm for solving the constrained form of Optimization Problem 10 .

3.2 Gradient Contraction for Constrained Optimization↩︎

Now let us consider the constrained form of the optimization problem given in 10 . To account for the constraint, we use a two-step gradient contractive algorithm which includes both a gradient step on variables \(x,\theta\) and a contraction step on variable \(\mathbf{u}\). In a standard gradient projection algorithm, we would evaluate the gradient with respect to \(x,\theta ,\mathbf{u}\) and then project \(\mathbf{u}\) onto the feasible set. However, because \(\mathbf{u}\) is infinite-dimensional and the feasible set is a manifold on this infinite-dimensional space, such an approach is difficult. As an alternative, however, one might propose Alg. [alg2:picard], which does not explicitly parameterize the variable \(\mathbf{u}\), but instead takes the gradient with respect to \(x\) and \(\theta\) for some fixed \(\mathbf{u}\) as \[\begin{align} \theta_{k+1} &= \Pi_\Theta[\theta_k - {\alpha \nabla_\theta L_{ls}(x_k, \theta_k, {{\mathbf{u}}}_k)}] \\ x_{k+1} &= \Pi_X[x_k - {\alpha \nabla_{x} L_{ls}(x_k, \theta_k, {{\mathbf{u}}}_k)}], \\[-8mm] \end{align}\] where note that, unlike in Alg. 1, there is no need to compute the gradient of a Picard iteration. Now, instead of parameterizing \(\mathbf{u}\) explicitly, for given \(x,\theta\), we use the Picard iteration to update our variable \(\mathbf{u}\) as \[{{{\mathbf{u}}}_{k+1}} = (1- \sigma) {{\mathbf{u}}}_k + \sigma \mathcal{P}_{ x_{k+1}, \theta_{k+1}} \mathbf{u}_k\] for some step size, \(\sigma \in (0, 1]\). Then, because the Picard iteration is a contraction on sufficiently short time intervals, this step not only updates \(\mathbf{u}\), but moves it closer to the feasible set \(\{\mathbf{u}\;:\; \mathbf{u} = {\mathcal{P}}_{x_{k+1}, \theta_{k+1}} \mathbf{u}\}\).

Now, Lipschitz continuity and strong convexity of the loss function with respect to \(x,\theta\) implies convergence of Alg. [alg2:picard] to a fixed point \(x^*_0\), \(\theta^*_0\) and \(\mathbf{u}^*_0\). However, our failure to explicitly parameterize \(\mathbf{u}\) does have a cost. Specifically, while the fixed point will be feasible, this fixed point will not necessarily be optimal. To see this, consider the KKT conditions necessary for optimality of Problem 10 .

Prop. [prop:optimality95conditions] implies the solution, \(x^*, \theta^*, \mathbf{u}^*\) of Problem 7 need satisfy Eqn. 12 . However, (interior) fixed points of Alg. [alg2:picard] only satisfy \(\nabla_{x, \theta} {L_{ls}}(x, \theta, \mathbf{u})=0\). To resolve this issue, in the following subsection, we approximate the multiplier \(\nabla_{x, \theta} \phi(t_i, x, \theta)\) in Eqn. 12 using Picard iteration and include this term in the updates of \(x\) and \(\theta\).

3.3 A Gradient Contraction Multiplier Algorithm for the Constrained Formulation↩︎

We now present Alg. [alg3:modified95picard], which does not require the use of infinite-dimensional variables and constraints, but for which the fixed point satisfies an approximation of the KKT conditions for optimality defined in Eqn. 12 of Prop. [prop:optimality95conditions]. Specifically, we require a fixed point \(\{x,\theta ,\mathbf{u}\}\) of Alg. [alg3:modified95picard] to satisfy \[\begin{align} &\nabla_{x,\theta} L_{ls}(x, \theta, \mathcal{P}^n_{ x, \theta}\mathbf{u}) = \nabla_{x,\theta} L_{ls}(x, \theta, \mathbf{u}) \label{eqn:KKT95approximation}\\[-1mm] &\quad + \sum\nolimits_{i = 1}^{N_s}\nabla_{u(t_i)} L_{ls}(x, \theta, \mathbf{u}) \nabla_{x,\theta} (\mathcal{P}^n_{ x, \theta}\mathbf{u})(t_i, x, \theta)=0\notag \\[-8mm]\notag \end{align}\tag{13}\] where recall \(\lim_{n\rightarrow \infty} \nabla_{x,\theta} ({\mathcal{P}}_{ x, \theta}^n \mathbf{u})(t) = \nabla_{x, \theta} \phi(t, x, \theta)\) as in Thm. [thm:picard95sensitivity]. Thus for any \(n\), the optimality conditions of Prop. [prop:optimality95conditions] are satisfied in some approximate sense, where the accuracy of this approximation increases with \(n\). As will be shown in Sec. 4 this allows us to prove convergence of the algorithm to optimality under suitable regularity and convexity conditions. As a practical matter, as will be shown in the numerical examples (Sec. 7), Alg. [alg3:modified95picard] requires only \(n=1\) or \(n=2\) for highly accurate solutions.

3.4 An Extended Gradient Contraction Multiplier Algorithm for the Constrained Formulation↩︎

Alg. [alg2:picard] and Alg. [alg3:modified95picard] eliminated the need for explicit parametrization of the variable \(\mathbf{u}\) by initializing \(\mathbf{u}\) with some arbitrarily chosen initial guess and then using the Picard operator in lieu of a projection to enforce the equality constraint \(\mathcal{P}_{x, \theta} \mathbf{u} - \mathbf{u}=0\). However, this approach is premised on the assumption that the Picard operator is sufficiently contractive so that, after a sufficient number of iterations, this equality constraint will be satisfied. However, the time interval for which the Picard iteration is contractive is substantially smaller than the interval on which the data is typically sampled or for which the solution map is defined. As a final step, therefore, we partition the time-domain into disjoint intervals and apply the Picard iteration to each interval, using a regularization term to ensure that discontinuities between intervals are minimized.

Specifically, for a given fixed interval of convergence, \(t \in [0,T]\), we lift the space of initial conditions as \(x \in X^{\otimes J}:=X \times \cdots \times X\) and the space of solutions as \(\mathbf{u} \in C[0,T]^{\otimes J}\) so that \(\mathbf{u}_{j}(0) = x_{j}\) and \(\mathbf{u}_{j}(T) = x_{j+1}\) for all \(j \in \overline{1, J-1}\) and \(\mathbf{u}_j(t)\) represents the estimated solution on interval \(t\in[(j-1)T, jT)\). Then the constrained optimization problem (Eqn. 7 ) may now be formulated as \[\begin{align} \min_{\substack{x \in X^{\otimes J},\, \theta \in \Theta\\ \mathbf{u} \in C[0,T]^{\otimes J}}} &\; \frac{1}{2N_s}\sum\nolimits_{i=1}^{N_s}\|y_i - g(t_i, \mathbf{u}_{\lfloor t_i/T \rfloor + 1}(t_i\bmod{T}), \theta)\|_2^2\notag\\[-2mm] \text{s.t.} \; & \mathcal{P}_{(j-1)T, x_j, \theta} \mathbf{u}_j = \mathbf{u}_j \qquad \forall j ,\, t\in [0,T] \label{eqn:optimization95problem95constrained95extended1} \\[-1mm] & \mathbf{u}_{j}(T) = x_{j+1} \qquad \forall j \in \overline{1, J-1}.\notag \\[-8mm] \notag \end{align}\tag{14}\] To see that Problems 10 and 14 are equivalent, we have the following.

To solve Optimization Problem 14 , we again use a Picard contractive step to avoid explicit parametrization of the variables \(\mathbf{u}_i\). However, this approach makes it difficult to directly enforce the linking constraints \(\mathbf{u}_{j}(T) = x_{j+1}\). Our approach, then is to move this constraint into the objective by using a penalty function as \[\begin{align} \label{eqn:optimization95problem95constrained95extended}\min_{\substack{x \in X^{\otimes J}, \theta \in \Theta\\ \mathbf{u} \in C[0,T]^{\otimes J}}} &\; \frac{1}{2N_s}\sum\nolimits_{i=1}^{N_s}\|y_i - g(t_i, \mathbf{u}_{\lfloor t_i/T \rfloor + 1}(t_i\bmod{T}), \theta)\|_2^2 \notag \\[-5mm] &\qquad\qquad + \lambda \sum\nolimits_{j=1}^{J-1}\|\mathbf{u}_{j}(T) - x_{j+1}\|_2^2\\ \text{s.t.} \;&\; \mathcal{P}_{(j-1)T, x_j, \theta} \mathbf{u}_j = \mathbf{u}_j \qquad \forall j \in \overline{1, J}, \notag\\[-8mm]\notag \end{align}\tag{15}\] where \(\lambda \geq 0\) is a regularization parameter.

Having formulated the constrained optimization problem, we now propose an extended version of the algorithm discussed in Subsec. 3.3 where Picard iterations are introduced in the objective function so that the gradient of this objective approximates the KKT conditions in Eqn. 12 – similar to the approach described in Eqn. 13 . Specifically, we now have objective function \[\begin{align} \label{eqn:extended95loss} & L_{\lambda, n}(x, \theta, \mathbf{u}) :=\sum_{j=1}^{J-1} \lambda \|\mathcal{P}^n_{(j-1)T, x_j, \theta} \mathbf{u}_{j}(T) - x_{j+1}\|_2^2\\[-3mm] & +\frac{1}{2N_s} \sum_{i=1}^{N_s}\|y_i - g(t_i, (\mathcal{P}^n_{(j-1)T, x_j, \theta} \mathbf{u}_{\lfloor t_i/T \rfloor+ 1})(t_i\bmod{T}), \theta)\|_2^2 \notag \\[-8mm] \notag \end{align}\tag{16}\] for which we define Alg. [alg4:ext95picard].

Although this approach allows us to approximate the solution of a parameter estimation problem on an arbitrary time interval, there are some technical challenges. First, the approach requires additional variables \(x_{j}\) (initial conditions at time points \((j-1)T\)) – increasing the dimension of the optimization problem. Second, the use of a penalty function in lieu of the constraint \(\mathbf{u}_j(T) = x_{j+1}\) implies discontinuities unless \(\lambda\) is chosen to be very large. However, increasing the value of \(\lambda\) increases the gradient of the objective – resulting in the need for smaller step sizes and hence slowing the convergence rate of the algorithm. These numerical issues are analyzed in detail in Subsec. 7.1.

4 Gradient-Contractive Algorithms↩︎

In this section, we propose sufficient conditions for convergence of the gradient-contractive algorithms presented in Sec. 3. Specifically, we propose a general class of optimization problems which includes those posed in 10 and 15 . Then, based on this generalized form of optimization problem, we define a generalized class of gradient-contractive algorithms. These two-step algorithms repeatedly update the optimization variables using both a gradient descent step and a contraction map (Algs. [alg2:picard], [alg3:modified95picard] and [alg4:ext95picard] are special cases of this approach). Finally, in Thm. [thm:gradient-contractive-map-is-contraction], we provide sufficient conditions such that each step of the gradient-contractive algorithm is a contraction – implying convergence.

To begin, suppose we are given a general optimization problem of the form \[\label{eqn:general95optimization95problem} \min_{\substack{x \in X, \mathbf{u} \in U}} \; L(x, \mathbf{u})\quad \text{s.t.} \;\; \mathcal{P}(x, \mathbf{u}) = \mathbf{u}, \tag{17}\] where \(U \subset \mathcal{F}(\mathbb{R}^{n_i}, \mathbb{R}^{n_u})\), \(X \subset \mathbb{R}^{n_x}\), \(\mathcal{P}:X\times U \rightarrow U\) and \(L:X\times U \rightarrow \mathbb{R}\). Now the gradient-contractive algorithm for step sizes \(\alpha > 0\), \(\sigma \in (0, 1]\) and order \(n\in \mathbb{N}\) is defined by the sequence \((x_{k},\mathbf{u}_{k})=\mathcal{T}^k (x_{0},\mathbf{u}_{0})\) where we say \((x_{k+1},\mathbf{u}_{k+1})=\mathcal{T} (x_{k},\mathbf{u}_{k})\) if \[\begin{align} \label{eqn:contractive95map} x_{k+1} &= \Pi_X\big[x_k - \alpha \nabla_x L(x_k, \mathcal{P}^n(x_k, \mathbf{u}_k))\big] \notag\\ \mathbf{u}_{k+1} &= (1 - \sigma) \mathbf{u}_{k} + \sigma \mathcal{P}(x_{k+1}, \mathbf{u}_k), \\[-8mm] \notag \end{align}\tag{18}\] where recall \(\Pi_X\) is the projection operator.

Clearly, Algs. [alg3:modified95picard] and [alg4:ext95picard] are special cases of this generalized gradient contractive algorithm, where \(n=0\) for Alg. [alg2:picard], \(\mathcal{P}\) is the Picard operator, \(L_i\) are as defined in Eqn. 6 and parameter variables \(\theta\) are combined into \(x\). Now, we provide conditions under which \(\mathcal{T}\) is a contraction, defined as follows.

When proving contractivity, we use the norm \(\lVert{x,\mathbf{u}}\rVert_H=\lVert{x}\rVert_X+\lVert{\mathbf{u}}\rVert_U\).

The key property of contractive mappings is that, when iteratively applied to any \(u\in H\), they converge to a fixed point, \(u^*=\lim_{k \rightarrow \infty} \mathcal{T}^k u_0\), for which \(\mathcal{T} u^*=u^*\).

Thm. [thm:gradient-contractive-map-is-contraction] gives sufficient conditions, expressed in terms of the properties of \(\mathcal{P}\) and \(L\), under which \(\mathcal{T}\) is a contraction and, which, therefore, implies that \(\mathcal{T}^k\) converges to a solution of Problem 17 .

Remark 1. The gradient-contractive approach of Thm. [thm:gradient-contractive-map-is-contraction] may be extended to other gradient-based algorithms (defined by alternative maps \(\mathcal{D}_n\)) as long as \(\mathcal{D}_n\) is contractive in \(x\) and Lipschitz in \(\mathbf{u}\).

To apply the results of Thm. [thm:gradient-contractive-map-is-contraction] to the parameter estimation problem, we define the mapping \(\mathcal{P}\) as a Picard iteration. In the following section, we will recall properties of the Picard iteration. Then, in Sec. 6, we will give conditions on the parameterized vector field, \(f\), under which the proposed modified and extended gradient contractive algorithms converge.

5 Contraction and Lipschitz continuity of the Picard Operator↩︎

In Sec. 4, we proposed a generalized form of gradient-contractive algorithm (Eqn. 18 ) for solving the class of optimization problems defined in Eqn. 17 , which are defined in terms of a Lipschitz separable objective, \(L(x, \mathbf{u})\) and an equality constraint of the form \(\mathcal{P}(x,\mathbf{u})=\mathbf{u}\). Furthermore, in Thm. [thm:gradient-contractive-map-is-contraction], we provided conditions on \(\mathcal{P}\) and \(L\) under which the gradient-contractive algorithm converges. In this section, we return to the original parameter estimation problem, as formulated in Optimization Problems 10 and 15 , where \(\mathcal{P}\) is now the Picard operator. In this case, we show that Algs. [alg3:modified95picard] and [alg4:ext95picard] are special cases of the generalized gradient-contractive algorithm and that the Picard operator, \(\mathcal{P}_{t_0, x, \theta}\), satisfies the conditions of Thm. [thm:gradient-contractive-map-is-contraction]. Specifically, in Thm. [thm:picard95contraction] and Lem. [lem:picard95is95lipschitz95in95parameters] (Subsec. 5.1), we show that \(\mathcal{P}_{t_0, x, \theta}\) satisfies Condition 1) of Thm. [thm:gradient-contractive-map-is-contraction] and in Lem. [lem:picard95sensitivity95is95continuous] (Subsec. 5.2), we show that \(\mathcal{P}_{t_0, x, \theta}\) satisfies Condition 2). Next, in Sec. 6, we give conditions under which \(L_{\lambda,n}\) (as defined in Eqn. 16 ) likewise satisfies the conditions of Thm. [thm:gradient-contractive-map-is-contraction].

5.1 Picard Iterations as a Contraction Map↩︎

In this subsection, we examine the Picard operator and show it satisfies Condition 1) of Thm. [thm:gradient-contractive-map-is-contraction]. Specifically, recall from Definition [def:picard95operator] that for any \(t_0,x,\theta\) and parameterized vector field, \(f\), the Picard operator is defined as \[(\mathcal{P}_{t_0, x, \theta} \mathbf{u})(t) := x+ \int_{0}^{t} f(s + t_0, \mathbf{u}(s), \theta) ds.\]

The contractive properties of the Picard operator and its gradient are relatively well-established on a sufficiently short time interval. However, in our case, we require these contractive properties to be uniform on domains \(x\in X\) and \(\theta \in \Theta\) which results in a slight variation on the standard result.

Thm. [thm:picard95contraction] shows that \(\mathcal{P}_{t_0, x, \theta}\) satisfies Condition 1) of Thm. [thm:gradient-contractive-map-is-contraction]. Before moving on to prove that Condition 2) is satisfied (in Subsec. 5.2), we state two further results which will be useful in establishing this result.

First, we note that, trivially, the solution map, \(\phi\), itself is a fixed point of the Picard iteration, which is formally stated as follows.

Second, we establish a uniform Lipschitz continuity bound on \(\mathcal{P}_{t_0, x, \theta}\).

Using these results, we now show that Condition 2) of Thm. [thm:gradient-contractive-map-is-contraction] is satisfied.

5.2 Lipschitz continuity of Gradients of Picard Iterations↩︎

In this subsection, we show that \(\mathcal{P}\) satisfies Condition 2) of Thm. [thm:gradient-contractive-map-is-contraction]. First, we use an extended version of the result in [32] (Thm. 6.20) to show that if the parameterized vector field is sufficiently smooth, the Picard operator, \(\mathcal{P}\) is differentiable in \(x,\theta\).

Iteratively applying Lem. [lem:picard95is95continuous], we have that \(\mathcal{P}^n\) is continuously differentiable with respect to \(t_0,x,\theta \in \Gamma \times X\times \Theta\). If \(f \in C^1\), we may thus define the gradients of an \(n\)-th order Picard operator as follows. \[S_n(t, x, \theta) = \begin{bmatrix}\nabla_x ({\mathcal{P}}_{t_0, x, \theta}^n \mathbf{u})(t) & \nabla_\theta ({\mathcal{P}}_{t_0, x, \theta}^n \mathbf{u})(t)\end{bmatrix}. \label{eqn:picard95gradient} \tag{20}\]

Next, we show that \(\{S_n\}_{n=1}^{\infty}\) is a Cauchy sequence in \(C([0, T]\times X \times \Theta)\) and hence converges to the gradient of the solution map, \(\phi\).

An immediate consequence of Lem. [lem:picard95sensitivity95is95cauchy] is that \(\{S_n\}_{n=1}^{\infty}\) converges to the gradient of the solution map.

Finally, we show that \(\nabla_{x, \theta} \mathcal{P}_{t_0, x, \theta}\) is Lipschitz continuous, thereby satisfying Condition 2) of Thm. [thm:gradient-contractive-map-is-contraction] (wherein through some abuse of notation, \(x \in X\) becomes \((x,\theta) \in X \times \Theta\) and \(U\) becomes \(C_a[0, T]\)).

To conclude, we have shown that \(\mathcal{P}_{t_0, x, \theta}\) satisfies conditions 1) and 2) of Thm. [thm:gradient-contractive-map-is-contraction]. Next, in Sec. 6, we give conditions under which \(L_{\lambda,n}\) (as defined in Eqn. 16 ) likewise satisfies the conditions of Thm. [thm:gradient-contractive-map-is-contraction] and use Thm. [thm:gradient-contractive-map-is-contraction] to show that Algs. [alg3:modified95picard] and [alg4:ext95picard] converge.

6 Conditions for Convergence of Gradient-Contractive Algorithm↩︎

In this section, we give conditions under which Algs. [alg3:modified95picard] and [alg4:ext95picard] converge to a fixed point. Specifically, we consider Algs. [alg3:modified95picard] to be a special case of Alg. [alg4:ext95picard] and use the results of Sec. 5 to propose conditions on the loss function, \(L_{\lambda, n}\) (as defined in 16 ), under which the conditions of Thm. [thm:gradient-contractive-map-is-contraction] are satisfied.

First, recall that the gradient-contractive algorithm for step sizes \(\alpha > 0\), \(\sigma \in (0, 1]\) and order \(n\in \mathbb{N}\) is defined by the sequence \((x_{k},\mathbf{u}_{k})=\mathcal{T}^k (x_{0},\mathbf{u}_{0})\) where we say \((x_{k+1},\mathbf{u}_{k+1})=\mathcal{T} (x_{k},\mathbf{u}_{k})\) if \[\begin{align} x_{k+1} &= \Pi_X\big[x_k - \alpha \nabla_x L(x_k, \mathcal{P}^n(x_k, \mathbf{u}_k))\big] \\ \mathbf{u}_{k+1} &= (1 - \sigma) \mathbf{u}_{k} + \sigma \mathcal{P}(x_{k+1}, \mathbf{u}_k). \\[-8mm]\notag \end{align}\]

Note that Alg. [alg4:ext95picard] is a form of gradient-contractive algorithm, with step sizes \(\alpha > 0\), \(\sigma \in (0, 1]\) and order \(n\in \mathbb{N}\), wherein \(x\in X\) becomes \((x,\theta) \in X^{\otimes J} \times \Theta\) and \(\mathbf{u} \in U\) becomes \(\{\mathbf{u}_j\}_{j=1}^J \in C[0,T]^{\otimes J}\) and where we say \(([x_{k+1},\theta_{k+1}],\mathbf{u}_{k+1})=\hat{\mathcal{T}}([x_{k},\theta_k],\mathbf{u}_{k})\) if \[\begin{align} \label{eqn:picard-gradient-contractive-iterations} x_{k+1, j} &= \Pi_{X}[x_{k, j} - \alpha \nabla_{x_{k, j}} L_{\lambda, n}(x_{k}, \theta_{k}, \mathbf{u}_{k})] \notag\\ \theta_{k+1} &= \Pi_{\Theta}[\theta_{k} - \alpha \nabla_{\theta} L_{\lambda, n}(x_{k}, \theta_{k}, \mathbf{u}_{k})] \\ \mathbf{u}_{k+1, j} &= (1 - \sigma) \mathbf{u}_{k, j} + \sigma \mathcal{P}_{(j-1)T, x_{k+1, j}, \theta_{k+1}} \mathbf{u}_{k, j}\notag. \\[-8mm]\notag \end{align}\tag{21}\]

We now show that if \(f,g\) are smooth and \(L_{\lambda,n}\) are strongly convex, Alg. [alg4:ext95picard] converges to a fixed point for some \(\sigma, \alpha\).

Remark 1. Note that in practice we have four hyper-parameters \(\alpha, \sigma, T\) and \(n\). One option to choose these parameters is to fix \(\sigma = \min\{1, \frac{\mu}{K_L - \mu}\}\) and \(\alpha = \frac{2}{\mu + K_L}\), where \(\mu\) is a strong convexity modulus of \(L_{\lambda, n}\) and \(K_L\) is Lipschitz constant of \(\nabla_{x, \theta} L_{\lambda, n}\). Then, for a fixed \(T\), \(n\) may be chosen sufficiently large. Alternately, for a fixed \(n \geq 1\) we may choose \(T\) sufficiently small.

7 Numerical experiments.↩︎

Having defined a gradient contractive algorithm based on the Picard mapping (Alg. [alg4:ext95picard]), and having established conditions for convergence of this algorithm to a solution of the parameter estimation problem, we now proceed to evaluate practical aspects of the performance of this algorithm. This analysis consists of five tests, each focusing on one metric of performance of the algorithm, and each using a nonlinear model and a data set indicative of that metric of performance. For our numerical results, the metric of performance is the average relative error (ARE) in the parameter estimates. Specifically, given \(N_{exp}\) estimates of the parameter values, ARE is defined as \[\label{eqn:ARE} ARE = \frac{1}{N_{exp}}.\sum_{i = 1}^{N_{exp}}\frac{| \hat{\theta}_{i} - \theta|}{|\theta|},\tag{22}\] where \(\hat{\theta}_i\) are the estimated parameters for the \(i^{th}\) data set and \(\theta\) is the true parameter values.

In all examples, initial conditions are unknown, and measurements include varying levels of noise. All tests of Alg. [alg4:ext95picard] use step size \(\sigma=.5\) and a backtracking line search to determine step size \(\alpha\). The parameter \(T\) varies by example and termination occurs when the norm of gradient \(\lVert{\nabla_{\theta}L_{\lambda,n}}\rVert\) is sufficiently small. The numerical tests are summarized as follows.

First, in Subsec. 7.1, we evaluate the effect of the regularization parameter \(\lambda\) on error in the parameter estimates using various levels of noise and as applied to the Van der Pol Oscillator, with direct measurement of both states. Second, in Subsec. 7.2, we evaluate the effect of irregular sampling intervals on error in parameter estimates and as illustrated by the FitzHugh-Nagumo neuron using measurements of a single state. Third, in Subsec. 7.3, we consider performance when data is obtained using multiple instances but with sparse data (long periods between samples) as illustrated by the Rosenzweig-MacArthur model. Fourth, in Subsec. 7.4 we consider experiments driven by external excitation and with a large number of states and parameters. We evaluate the effect of different excitations on estimation accuracy as illustrated by a model of tumor growth. Finally in Subsec. 7.5, we compare the accuracy of the algorithm with both comparable heuristics [8] and gradient-free optimization [23], [33] as illustrated using the Lorenz model. Additional numerical examples for estimating unknown parameters in real-world models of the attitude dynamics can be found in [34].

7.1 The Van der Pol Oscillator and the Effect of Regularization Parameter↩︎

In the first numerical example, we investigate the performance of Alg. [alg4:ext95picard] (defined by regularized loss function \(L_{\lambda, n}\) in Eqn. 16 ) as a function of regularization parameter \(\lambda\). Recall, this regularization parameter \(\lambda\) weights a penalty for discontinuity of the approximated solution at times \(jT\). Hence, for a larger weight \(\lambda\), the estimated solution, \(\mathbf{u}\) is closer to the solution map on extended time intervals and hence better captures the relationship between estimated parameter values and evolution of the state. However, as \(\lambda\) increases, convergence time increases and less weight is placed on matching the estimated solution to the given data. To provide guidance on choice of regularization parameter, we consider Van der Pol oscillator. \[\begin{align} \dot{x}_1(t) &= \theta_1 x_2(t) \label{eqn:example95vdp95model} \\ \dot{x}_2(t) &= -\theta_1 x_1(t) + \theta_2 (1 - x_1^2(t)) x_2(t), \notag \\[-8mm]\notag \end{align}\tag{23}\] where both states are observable and our goal is to estimate the parameters \(\theta_1, \theta_2\). The data sets used in this example are generated from a single instance (with single initial condition \(x_0\)) over a time interval \([0,t_f]\) using \(t_f=10\) and true parameter values \(\theta_1 = \theta_2 = 1\), \(N_s\) evenly distributed sampling times \(t \in \{0:t_f/(N_s-1) i: t_f\}\) and noisy data of the form \(y_i=(1 + n_i)\phi(t_i, x_0, \theta) + m_i\), where \(\phi(t_i, x, \theta)\) are the actual states and \(n_i, m_i\) are normally distributed with zero mean and variance \(\varepsilon\).

To apply Alg. [alg4:ext95picard], we use the measurement data to approximate a bound on the Lipschitz constant as \(K_x\sim 4\) which yields a convergence interval of \(T = {1/K_x}= .25s\). This yields \(J={t_f/T}= 40\) sub-intervals. A single Picard iteration \(n=1\) is used.

To evaluate the effect of regularization, we now construct 100 data sets, \(Y_{1,\varepsilon,x_0}\) (using \(N_s=101\)) and \(Y_{2,\varepsilon,x_0}\) (using \(N_s=51\)) for variances \(\varepsilon \in \{0.01, 0.05, 0.1, 0.2, 0.4\}\) and 10 randomly selected initial conditions \(x_0\in X_{\lambda,\varepsilon}\) (such that \(x_0\) is always inside the limit cycle). Alg. [alg4:ext95picard] is then applied to each data set using regularization parameters \(\lambda \in \{1,2,5,10,15\}\) – yielding parameter estimates \(\hat{\theta}_{\lambda, \varepsilon,x_0}\) and initial condition estimates \(\hat{x}_{\lambda, \varepsilon,x_0}\).

For each choice of \(\lambda,\varepsilon\), performance is evaluated by first performing numerical simulation with estimated parameters, \(\hat{\theta}_{\lambda, \varepsilon,x_0}\), and initial conditions, \(\hat{x}_{\lambda, \varepsilon,x_0}\), and then computing the average loss so that if \(\phi\) is the solution map, we have \[L(\lambda, \varepsilon) = \frac{1}{10}\sum_{x_0\in X_{\lambda,\varepsilon}}\frac{1}{N_s} \sum\nolimits_{i=1}^{N_s}\|\phi(t_i,\hat{x}_{\lambda, \varepsilon,x_0},\hat{\theta}_{\lambda, \varepsilon,x_0}) - y_i\|_2^2.\] In addition, for a parameter, \(\theta\), we compute the average relative error (\(ARE\)) in the parameter estimates as in Eqn. 22 , where the summation is given by parameter estimates for different initial conditions \(x_0\in X_{\lambda,\varepsilon}\). As seen in Fig. 2, increasing the regularization parameter \(\lambda\) significantly improves both average loss and ARE for large noise variance (\(\varepsilon \in \{0.1, 0.2, 0.4\}\)), but does not significantly improve performance with low noise (\(\varepsilon \in \{0.01, 0.05\}\)). In addition, as expected, more data (\(N_s = 101\) vs \(N_s = 51\)) improved performance in almost every case – See Fig. 2 a, c) and b, d).

Although Alg. [alg4:ext95picard] accurately estimates parameters of the model in 23 , it is also possible to have locally identifiable or unidentifiable parameters. To illustrate the locally identifiable parameters, let’s consider a different parameterization of the model in 23 . Specifically, using a variable substitution \(\theta_1 = (\rho - 1)^2\), we have that \(\rho = 2\) and \(\rho = 0\) gives us \(\theta_1 = 1\), both choices leading to identical system response in 25 . Numerically, using Alg. [alg4:ext95picard] to estimate \(\rho\), we observe that initial parameter estimate \(\rho_{0}=-1\) results in convergence to \(\rho = 0\) and \(\rho_{0}=3\) results in convergence to \(\rho = 2\). This illustrates the importance of the choice of parameterization of the model to avoid locally identifiable or non-identifiable models [26].

7.2 The FitzHugh-Nagumo Neuron and the Effect of Irregular Sampling↩︎

In this subsection, we evaluate the effect of irregular sampling intervals on error in estimated parameters. For this case, we consider the FitzHugh-Nagumo model a single cortical neuron, which represents spiking and bursting behavior of neurons [35]. This model is a standard example of a stiff ODE and hence it is known that sampling at times of spiking and bursting are more useful than during interim periods – See Fig. 3. Specifically, we consider the model \[\begin{align} \dot{x}_1(t) &= 10(x_1(t) - 1/3 x_1^3(t) - x_2(t) + I(t)) \notag \\ \dot{x}_2(t) &= 10(\theta_1 x_1(t) - \theta_2 x_2(t) + 0.056) \label{eqn:FHN95example95model} \\ y(t) & = g(x(t)):=x_1(t), \notag \\[-8mm]\notag \end{align}\tag{24}\] where the states, \(x_1(t), x_2(t) \in \mathbb{R}\), represent membrane potential of neurons and membrane recovery variables and the input, \(I(t)\), represents the stimulus current and is set to \(I(t) = 1\) when the neuron is active. The measured output is \(y(t)\) and the parameters to be estimated are \(\theta_1, \theta_2 \geq 0\).

For this test, we use true parameter values \(\theta_1 = 0.08\) and \(\theta_2 = 0.064\). Based on observations of the output, \(y(t)=x_1(t)\), sampling times are clustered around observed spiking and bursting times, \(T_2:=\{1.55+0.0125k\}_{k=0}^{56}\) and \(T_4:= \{3.5+0.0125k\}_{k=0}^{40}\). Sampling times during recovery periods are less frequent: \(T_1:= \{0.25k\}_{k=0}^6\), \(T_3 := \{2.3+0.25k\}_{k=0}^4\), and \(T_5:= \{4.05 + 0.25k\}_{k=0}^3\) where the union of all sampling times is defined as \(t_i \in T := \cup_{i=1:5} T_i\). These sampling times are depicted as blue dots in Fig. 3 and where (as in Subsec. 7.1) measurements are taken over a single instance using initial condition \(x_1(0)=0, x_2(0)=2\) as \[y_i=(1 + n_i)g(\phi(t_i, x(0), \theta)) + m_i\qquad t_i \in T,\] where \(\phi\) is the solution map and \(n_i, m_i\) are normally distributed with zero mean and variances of \(\varepsilon = .05\) (Fig 3a) and \(\varepsilon = .2\) (Fig 3b). Applying Alg. [alg4:ext95picard] with \(\lambda =5\), \(n =2\), and \(T=.1\) we obtain parameter estimates \(\hat{\theta}_1 = 0.079\) and \(\hat{\theta}_2 = 0.063\) for variance \(\varepsilon = .05\); and \(\hat{\theta}_1 = 0.087\) and \(\hat{\theta}_2 = 0.067\) for variance \(\varepsilon=.2\).

7.3 Rosenzweig-MacArthur predator-prey model and the Effect of Sparse Data.↩︎

In the next example, we evaluate performance of Alg. [alg4:ext95picard] for sparse data using the Rosenzweig-MacArthur predator-prey model [36], wherein the states represent populations of predators and times are measured in years. As in many biological models, collecting data entails the labor-intensive counting of populations and hence data for such models is typically limited to a few samples. In addition, measurements of biological processes are often inaccurate and hence repeated experiments are used to obtain statistically significant sample sizes. For example, in the case of predator-prey models, data may be based on tracking several distinct populations in order to reduce the influence of unmodelled and localized factors such as geography. In our case, the model takes the following form [36] \[\begin{align} \dot{x}_1(t) &= rx_1(t)\big(1 - \frac{x_1(t)}{K}\big) -\frac{ ax_1(t)x_2(t)}{b + x_1(t)}\notag \\[-2mm]\dot{x}_2(t) &= k \frac{ ax_1(t)x_2(t)}{b + x_1(t)} - c x_2(t), \label{eqn:example95preypredator95model}\\[-9mm]\notag \end{align}\tag{25}\] where \(x_1(t), x_2(t)\) represent prey and predator populations, respectively. This model has \(6\) unknown parameters (\(a,b,c,r,k,K\)), with: \(r, K\) representing reproduction rate and carrying capacity of the prey population; \(a, k\) are encounter and growth rates, \(b\) is an environmental factor, and \(c\) is predator death rate. For this test, we use true parameters \(a = 2.8\), \(b=0.7\), \(c=1.35\), \(r = 3.5\), \(k = 1.5\) and \(K = 1.4\).

To evaluate the effect of sparse sampling and multiple instances, the algorithm is tested using data obtained from between 1 and 4 distinct populations (\(N_c=1:4\) is number of initial conditions), with population counts being made at years \(t_i \in \{0, 1.25, 2.5, 3.75, 5\}\) using measurement with noise model \(y_{ij} = (1 + n_i) \phi(t_i, x_j, \theta) + m_i\), where \(n_i, m_i\) are normally distributed with zero mean and variance \(\varepsilon = 0.05\).

Alg. [alg4:ext95picard] is then applied to each case with parameters \(\lambda = 5\), \(n = 2\) and \(T = 0.2\). The test is repeated ten times and the average relative error in the parameter estimates is listed in Table ¿tbl:tab:pp95results?, where even for small amount of samples we were able to identify parameters.

As expected, these results demonstrate that increasing the number of sampled populations reduces error in the parameter estimates, with reasonable accuracy achieved using as few as two populations. However, certain parameters (\(a\) and \(b\)) do not appear to converge, even with an increased number of populations. Furthermore, we find that for noise variance \(\varepsilon = 0.05\), the objective value for, e.g., \(a=3.65\) and \(b=0.83\) is indistinguishable from the objective value obtained using the true values of \(a=2.8\) and \(b=.7\). This implies that the problem here is not the existence of local minima, but rather that the model parameterization may not be identifiable as stated.

7.4 Tumor Growth, Identifiability, and Excitation.↩︎

In this example, we examine the significance of model excitation for a weakly identifiable model. As discussed in [37], [38], biological models with multiple states and parameters are often weakly identifiable in that multiple choices of parameters yield the same or approximately the same solution. We examine here the effect of using datasets generated from multiple non-redundant excitations to improve identifiability of these parameters. Specifically, we use a model of tumour immunodynamics with chemotherapy [39] where different drug dosing strategies are used as inputs to the model, given as \[\begin{align} \dot{N}(t) &= \theta_1 N(t)(1 -N(t)) -\theta_1 N(t) T(t) -(1-e^{-v(t)}) N(t) \notag\\ \dot{T}(t) &= \theta_2 T(t)(1 - T(t)) - \theta_3 T(t) I(t) - \theta_1 T(t) N(t) \notag\\[-2mm] & \qquad\qquad\qquad\qquad\qquad\qquad- 3(1 - e^{-v(t)}) T(t)\notag\\[-2mm] \dot{I}(t) &= \theta_4 + \theta_5 \frac{I(t) T(t)}{\rho + T(t)} - \theta_{1} I(t) T(t) - \theta_{6} I(t) \label{eqn:example95tumor95model}\\[-4mm] & \qquad\qquad\qquad\qquad\qquad\qquad - 2 (1 - e^{-v(t)}) I(t) \notag \\[-2mm] \dot{v}(t) &= - v(t) + u(t) \notag \\[-2mm] y(t) &= g(N(t), T(t), I(t), v(t)) := \begin{bmatrix}N(t) & T(t) & I(t)\end{bmatrix}^T, \notag\\[-8mm]\notag \end{align}\tag{26}\] where \(N(t)\) is the density of healthy cells, \(T(t)\) is tumor cells, \(I(t)\) is immune cells and \(v(t)\) is the concentration of drugs in the tumor area. The system has a single input \(u(t)\), that represents drug dosing. We assume that states \(N(t), T(t), I(t)\) are outputs of the system. Parameter values are taken from [39], where \(\theta_{1} = 1\), \(\theta_2 = 1.5\), \(\theta_3 = 0.5\), \(\theta_4 = 0.33\), \(\theta_5 = 0.01\), \(\theta_{6} = 0.2\) and \(\rho = 0.3\). The initial conditions are \(N(0) = 1\), \(T(0) = 2\), \(I(0) = 1.65\), \(v(0) = 0\).

5 drug dosing strategies are chosen as \(u_1(t) = 0\), \[\begin{align} u_2(t) & = \begin{cases} e^{-t}, & \text{if } t \in [0, 3] \\ 0, & \text{otherwise} \end{cases} \quadu_3(t) = \begin{cases} e^{t}, & \text{if } t \in [0, 3] \\ 0, & \text{otherwise} \end{cases}\\ u_4(t) & = \begin{cases} 3 - t & \text{if } t \in [3, 6] \\ 0, & \text{otherwise} \end{cases} \quadu_5(t) = \begin{cases} e^t, & \text{if } t \in [1.5, 4.5] \\ 0, & \text{otherwise,} \end{cases} \\[-9mm]\notag \end{align}\] where it is observed that \(u_3,u_4,u_5\) clear the tumour, but \(u_1\) and \(u_2\) do not.

For each dosing strategy, \(i=1,\cdots,5\), we generate 3 datasets \(\{Y_{i,j}\}_{j=1}^3\) each using measurement noises \(n_i, m_i\) (normally distributed with zero mean and variance \(\varepsilon=.01\)) so that \(y_{i} = (1 + n_i)g(\phi_u(t_i, x, \theta)) + m_i\) where \(\phi_u\) is the solution of Eqn. 26 for input \(u(t)\), and \(t_i = \{0:\frac{t_f}{N_s - 1}i : t_f\}\), where \(N_s = 51\) and \(t_f = 6\).

For testing, we apply Alg. [alg4:ext95picard] (with \(\lambda = 2\), \(n = 2\) and \(T = 0.1\)) to the datasets \(Y_{1,j}\), \(Y_{1,j}\cup Y_{2,j}\), \(Y_{1,j}\cup Y_{2,j}\cup Y_{3,j}\), \(Y_{1,j}\cup Y_{2,j}\cup Y_{3,j}\cup Y_{4,j}\), and \(Y_{1,j}\cup Y_{2,j}\cup Y_{3,j}\cup Y_{3,j}\cup Y_{5,j}\) – representing an increase in the number of excitations used to identify the data. In each case we determine the ARE (over \(j=1:3\)) for each identified parameter. The results are listed in Table 1, where green indicates that the ARE for the given parameter is less that \(0.3\). These results indicate the failure to include any excitation (\(Y_{1j}\)) results in failure to effectively identify the parameters while inclusion of inputs \(u_2,u_3,u_4\) generates improved results. Inclusion of all 5 excitations generates the best results, although even in this case, it seems that some parameters are not identifiable. This result suggests that the algorithm converges to local optima of the loss function, a typical behavior for parameter estimation problems. Indeed, the achieved objective value for the estimated parameters is \(L_{ls} = 0.027\) compared to \(L_{ls} = 0.006\) for the true parameter values. Thus, in contrast to Example 7.3, the proposed algorithm was not able to find a globally-optimal solution. The failure of the convergence may be caused by the non-convexity of the optimization problem or the convergence to approximated KKT conditions. Thus, to more accurately identify the parameters, reparameterization of the model may be necessary (Sec. 5 in [40] or Sec. 2 in [37]).

7.5 The Lorenz System: Chaos and Comparison.↩︎

The goal of this last numerical test is to evaluate the accuracy of estimated parameters in the proposed algorithms as compared to existing state-of-the-art methods for parameter estimation. For this analysis we use data generated from the well-studied Lorenz system [2], [4], originally proposed as a model of convection rolls in the atmosphere [41]. Solutions to the Lorenz system tend to converge to a 2D manifold known as a “strange attractor”, but are chaotic within this manifold – implying adjacent trajectories diverge exponentially fast. As a result, while parameter estimation for the Lorenz system is reasonably accurate when measurement noise is small, many methods for parameter estimation fail when noise is significant. The Lorenz system is defined as \[\begin{align} \dot{x}_1(t) &= 10\sigma(x_2(t) - x_1(t)) \notag \\ \dot{x}_2(t) &= 10x_2(t)(\rho - x_3(t)) - x_2(t) \label{eqn:example95lorenz95model} \\ \dot{x}_3(t) &= 10x_1(t) x_2(t) - \beta x_3(t), \notag \\[-8mm]\notag \end{align}\tag{27}\] where \(x_1, x_2, x_3\) are states and \(\sigma, \rho, \beta\) are parameters.

For this analysis, we use Lorenz parameters \(\sigma = 1\), \(\beta = \frac{8}{3}\) and \(\rho = 2.8\). In order to allow for comparison with existing methods, each dataset uses a single initial condition, and directly measures the state at frequent, regular sampling times \(t_i = \{0: 2i/(N_s - 1) : 2\}\) (\(N_s = 101\)). Measurement noise is modelled as \(y_i = (1 + n_i)\phi(t_i, x_0, \theta) + m_i,\) where \(n_i, m_i\) are normally distributed with zero mean and variance \(\varepsilon\). This approach is used to obtain \(50\) data sets \(Y_{\varepsilon, x_0} = \{(y_i, t_i) \}_{i = 1}^{N_s}\), where \(\varepsilon \in \{0, 0.01, 0.05, 0.1, 0.2, 0.4\}\) and 10 randomly chosen initial conditions \(x_0\), normally distributed with zero mean and variance \(1\).

Parameter estimation is then performed on each dataset using each of 6 algorithms: Picard (Alg. [alg4:ext95picard]); SINDy [2]; fmincon [16]; GA [23]; AMIGO2 [42]; and SciML Sensitivity Toolbox [43]. For Picard, Alg. [alg4:ext95picard] is used with parameters \(\lambda = 5\), \(n = 2\) and \(T = 0.05\). For SINDy, we use a slight modification of the black-box implementation given in [2] (regularization parameter \(\lambda=0.25\) and measurements of \(\dot{x}(t_i)\) are computed using first-order difference). For fmincon, we use nlgreyest from MATLAB’s System Identification Toolbox [33] with selected optimizer: fmincon. For AMIGO2 Toolbox, we use Enhanced Scatter Search to find the optimal parameters [44]. For SciML Sensitivity Toolbox, we use Adam optimizer followed by BFGS optimizer [45]. For GA, we designed a gradient-free black-box optimization approach to solution of Eqn. 6 . Specifically, for given \(x,\theta\), we use numerical simulation to compute the solution \(\phi\) which then can be used to compute the least squares loss. With this loss oracle in hand, we use a genetic algorithm (MATLAB’s ga function – with tolerances \(0.01\) and run time of \(5\) mins on Intel i7-4960X CPU at 3.60 GHz) to solve the optimization problem. For all methods, parameters are restricted to \(\rho \in [1.4 , 4.2]\), \(\sigma \in [0.5 ,1.5]\), \(\beta \in [\frac{4}{3} , 4]\), and \(x_i(0)=[0.5 x_i , 1.5 x_i]\), where \(x_i\) are true initial conditions.

Each of the four algorithms was evaluated on each of the 50 datasets and for each \(\varepsilon\), we computed ARE and standard deviation as described previously and illustrated in Fig 4. These results show that the Picard method outperforms all methods for low level of noise \(\varepsilon \leq 0.1\) and for larger values of noise (\(\varepsilon \geq 0.2\)) Picard performs significantly better than SINDy and fmincon methods – but comparable to the proposed black box optimization alternative to the gradient-based Picard approach. Simulation of the Lorenz system and comparison to data for estimated parameters in each method (with noise variance \(\varepsilon=0.4\) and single dataset) is given in Fig. 5.

8 Conclusion↩︎

In this paper, we reformulate the problem of parameter estimation for nonlinear ordinary differential equations as a constrained optimization problem with infinite-dimensional variables and constraints. We then propose a new class of gradient-contractive algorithms, based on the contractive properties of the Picard iteration, which eliminates the need for infinite-dimensional constraints and variables (Algs. [alg2:picard]-[alg4:ext95picard]). In contrast with existing methods, this approach is gradient-based and allows for data which is sparse, irregular, and includes only partial measurements of the state and does not require numerical simulation or measurements of time derivatives. Furthermore, we have proposed sufficient conditions under which the proposed algorithm converges to a local optima. Finally, we have exhaustively tested several aspects of performance of the algorithm on a battery of nonlinear models and noisy datasets.

Acknowledgments↩︎

We would like to thank Rolf Findeisen for his feedback and acknowledge the priority of his team in recognizing the role of Picard iteration in parameter estimation.

9 Auxiliary Lemmas↩︎

In this appendix, we prove the auxiliary lemmas.

10 Gradients of Picard Iterations↩︎

In this appendix, we show that gradients of Picard Iterations are Lipschitz continuous. Specifically, we first show that \(\nabla_{x, \theta} \mathcal{P}_{t_0, x, \theta} \mathbf{u}\) is Lipschitz continuous with respect to \(\mathbf{u}\) uniformly in \((x, \theta) \in X\times \Theta\).

Next, we show that \(\nabla_{x, \theta} \mathcal{P}_{t_0, x, \theta} \mathbf{u}\) is Lipschitz continuous with respect to \(x\) and \(\theta\) uniformly in \(\mathbf{u} \in C_a[0, T]\).

References↩︎

---

1. This work was supported in part by the National Science Foundation under Grants NSF CCF-2323532 and NSF/NIH 2054354 awarded to M. M. Peet. (Corresponding author: Aleksandr Talitckii.)↩︎