Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Picard Iteration for Parameter Estimation in Nonlinear Ordinary Differential Equations

Aleksandr Talitckii, Matthew M. Peet

TL;DR

This work develops a Picard-operator–based constrained reformulation for parameter estimation in nonlinear ODEs when data are noisy, sparse, and irregular, eliminating the need for an explicit solution map. It introduces gradient-contractive algorithms that couple gradient steps on the parameters with contraction steps on the infinite-dimensional state trajectory, underpinned by contraction and Lipschitz properties of Picard iterations. The authors prove convergence to a local optimum under suitable regularity and convexity, and validate the approach on multiple nonlinear models, including Van der Pol, FitzHugh–Nagumo, Rosenzweig–MacArthur, tumor growth, and Lorenz systems, highlighting robustness to irregular sampling, sparse data, and measurement noise. The framework supports long horizons via extended Picard iterations and interval-based decomposition, with practical guidance on parameter choices and regularization. Overall, the method offers a principled, data-efficient alternative to black-box optimization for identifying nonlinear dynamics from imperfect observations, with demonstrated applicability to identifiability and excitation design problems.

Abstract

We consider the problem of using experimental time-series data for parameter estimation in nonlinear ordinary differential equations, focusing on the case where the data is noisy, sparse, irregularly sampled, includes multiple experiments, and does not directly measure the system state or its time-derivative. To account for such low-quality data, we propose a new framework for gradient-based parameter estimation which uses the Picard operator to reformulate the problem as constrained optimization with infinite-dimensional variables and constraints. We then use the contractive properties of the Picard operator to propose a class of gradient-contractive algorithms and provide conditions under which such algorithms are guaranteed to converge to a local optima. The algorithms are then tested on a battery of models and variety of datasets in order to demonstrate robustness and improvement over alternative approaches.

Picard Iteration for Parameter Estimation in Nonlinear Ordinary Differential Equations

TL;DR

This work develops a Picard-operator–based constrained reformulation for parameter estimation in nonlinear ODEs when data are noisy, sparse, and irregular, eliminating the need for an explicit solution map. It introduces gradient-contractive algorithms that couple gradient steps on the parameters with contraction steps on the infinite-dimensional state trajectory, underpinned by contraction and Lipschitz properties of Picard iterations. The authors prove convergence to a local optimum under suitable regularity and convexity, and validate the approach on multiple nonlinear models, including Van der Pol, FitzHugh–Nagumo, Rosenzweig–MacArthur, tumor growth, and Lorenz systems, highlighting robustness to irregular sampling, sparse data, and measurement noise. The framework supports long horizons via extended Picard iterations and interval-based decomposition, with practical guidance on parameter choices and regularization. Overall, the method offers a principled, data-efficient alternative to black-box optimization for identifying nonlinear dynamics from imperfect observations, with demonstrated applicability to identifiability and excitation design problems.

Abstract

We consider the problem of using experimental time-series data for parameter estimation in nonlinear ordinary differential equations, focusing on the case where the data is noisy, sparse, irregularly sampled, includes multiple experiments, and does not directly measure the system state or its time-derivative. To account for such low-quality data, we propose a new framework for gradient-based parameter estimation which uses the Picard operator to reformulate the problem as constrained optimization with infinite-dimensional variables and constraints. We then use the contractive properties of the Picard operator to propose a class of gradient-contractive algorithms and provide conditions under which such algorithms are guaranteed to converge to a local optima. The algorithms are then tested on a battery of models and variety of datasets in order to demonstrate robustness and improvement over alternative approaches.
Paper Structure (21 sections, 89 equations, 5 figures, 2 tables, 4 algorithms)

This paper contains 21 sections, 89 equations, 5 figures, 2 tables, 4 algorithms.

Table of Contents

  1. Introduction
  2. Formulating the Parameter Estimation Problem for Nonlinear Differential Equations
  3. Reformulation using the Picard Operator
  4. Gradient Descent for Parameter Estimation
  5. A Gradient Descent Algorithm for the Unconstrained Formulation
  6. Gradient Contraction for Constrained Optimization
  7. A Gradient Contraction Multiplier Algorithm for the Constrained Formulation
  8. An Extended Gradient Contraction Multiplier Algorithm for the Constrained Formulation
  9. Gradient-Contractive Algorithms
  10. Contraction and Lipschitz continuity of the Picard Operator
  11. Picard Iterations as a Contraction Map
  12. Lipschitz continuity of Gradients of Picard Iterations
  13. Conditions for Convergence of Gradient-Contractive Algorithm
  14. Numerical experiments.
  15. The Van der Pol Oscillator and the Effect of Regularization Parameter
  16. ...and 6 more sections

Figures (5)

  • Figure 1: Loss function as in Eqn. \ref{['eqn:optimization_problem_unconstrained']} and normalized error of parameters ($ARE$) for Van der Pol Oscillator in Eqn. \ref{['eqn:example_vdp_model']} identified by Alg. \ref{['alg4:ext_picard']} for $n = 1$, $T = 0.25$ and $\lambda \in \{1, 2, 5, 10, 15\}$. The true parameters, $\theta_1 = \theta_2 = 1$, variance of the measurements is $\varepsilon \in \{0.01, 0.05, 0.1, 0.2, 0.4\}$. The number of measurements is $N_s = 101$ or $N_s = 51$. The data for $N_s = 51$ is a subsample of the measurements $N_s = 101$ for twice larger sampling times. Average and standard deviation of $ARE$ of parameters and loss function is based on $10$ trials with randomized initial conditions.
  • Figure 2: Parameter Estimation using Alg. \ref{['alg4:ext_picard']} for the FitzHugh-Nagumo model in Eqn. \ref{['eqn:FHN_example_model']} with irregular sampling times (blue dots) clustered around spiking and bursting with noise variances of $\varepsilon=.05$ (Fig 2a) and $\varepsilon=.2$ (Fig 2b). The red dashed line indicates simulation based on the estimated parameters ($\hat{\theta}_i$) and initial states. True parameter values are $\theta_1 =0.08$ and $\theta_2 = 0.064$.
  • Figure 3: Accuracy in identified system parameters for the Lorenz system (Eqn. \ref{['eqn:example_lorenz_model']}) using Alg. \ref{['alg4:ext_picard']} (Picard), as compared with: SINDy; MATLAB's nlgreyest (w. fmincon) (fmincon); and a black box gradient-free alternative approach to solving Eqn. \ref{['eqn:optimization_problem_unconstrained']} (GA). For each algorithm, ARE and standard deviations for each parameter $\rho$ (Fig. 3a), $\sigma$ (Fig. 3b), $\beta$ (Fig. 3c) are computed over 10 datasets for each level of measurement noise variance -- $\varepsilon \in \{0.01, 0.05, 0.1, 0.2, 0.4\}$. True parameter values are $\rho = 2.8$, $\sigma = 1$, $\beta = 8/3$.
  • Figure 4: Numerical simulation of 3 states of the Lorenz system (Eqn. \ref{['eqn:example_lorenz_model']}) using identified parameters from Picard (Alg. \ref{['alg4:ext_picard']}), SINDy, fmincon and GA alternatives (See Fig. \ref{['fig:lorenz_results']} and description) for noise variance $\varepsilon=0.4$. True parameter values are $\rho =2.8$, $\sigma = 1$, $\beta = 8/3$; initial states are $x(0)\space =\space [1.46, 0.90, -1.57]^T$ ; and blue dots are measurements.
  • Figure : Gradient Descent Algorithm