# Multi-objective SINDy for parameterized model discovery from single
transient trajectory data

Javier A. Lemus  Benjamin Herrmann<EMAIL_ADDRESS>Department of
Mechanical Engineering, University of Chile, Beauchef 851, Santiago, Chile

###### Abstract

The sparse identification of nonlinear dynamics (SINDy) [1] has been
established as an effective technique to produce interpretable models of
dynamical systems from time-resolved state data via sparse regression.
However, to model parameterized systems, SINDy requires data from transient
trajectories for various parameter values over the range of interest, which
are typically difficult to acquire experimentally. In this work, we extend
SINDy to be able to leverage data on fixed points and/or limit cycles to
reduce the number of transient trajectories needed for successful system
identification. To achieve this, we incorporate the data on these attractors
at various parameter values as constraints in the optimization problem. First,
we show that enforcing these as hard constraints leads to an ill-conditioned
regression problem due to the large number of constraints. Instead, we
implement soft constraints by modifying the cost function to be minimized.
This leads to the formulation of a multi-objective sparse regression problem
where we simultaneously seek to minimize the error of the fit to the
transients trajectories and to the data on attractors, while penalizing the
number of terms in the model. Our extension, demonstrated on several numerical
examples, is more robust to noisy measurements and requires substantially less
training data than the original SINDy method to correctly identify a
parameterized dynamical system.

##  I Introduction

The behavior of complex systems found across science and engineering often
depends on quantities such as material or media properties, geometrical
parameters, the magnitude of the physical scales involved, boundary
conditions, and/or control inputs. The effect of parameter changes in
nonlinear dynamics has been the focus of intense research for decades due to
the emergence of bifurcations and pattern formation [2]. Moreover, reduced-
order models that capture the parameterized dynamics of complex engineering
systems [3, 4] are a critical enabler for real-time applications [5, 6, 7],
many-query optimization tasks [8, 9, 10, 11], and the development of digital
twins [12, 13, 14], which will have a profound impact on industries such as
aerospace [15, 16], electromobility [17], and biomedical [18].

Traditionally, the governing equations of dynamical systems have been derived
from first principles. However, although we are be able to build physics-based
models for many complex systems, these models may be computationally expensive
to simulate or may show poor predictive performance due to simplified physics
assumptions. Furthermore, governing equations derived from first principles
are not available in fields such as neuroscience, epidemiology, and finance.
Fortunately, the unprecedented availability and quality of data has led to the
emergence of data-driven model discovery — a paradigm in which we learn
governing laws directly from measurements of the system of interest [19]. In
recent years, a plethora of approaches have been proposed in response to this
trend, such as combined machine learning and data assimilation frameworks [20,
21], deep learning approaches [22, 23, 24, 25, 26, 27, 28], divide-and-conquer
strategies [29], symbolic regression via genetic algorithms [30, 31, 32],
Gaussian processes [33], nonlinear autoregressive algorithms [34], the dynamic
mode decomposition (DMD) [35, 36, 37] and its variants [38, 39, 40, 41],
Koopman theory [42, 43, 44], operator inference [45, 46], regression onto
spectral submanifolds [47, 48, 49], and the sparse identification of nonlinear
dynamics (SINDy) [1, 50]. The latter method, which is the focus of the present
work, has gained a considerable amount of traction because it follows a simple
philosophy, is easy to implement, and is highly extensible [1].

SINDy builds parsimonious and interpretable models of dynamical systems from
data by leveraging sparsity promoting regression techniques [51] that balance
model complexity and accuracy to avoid overfitting [1]. Since its
introduction, the method has been rapidly adopted to address challenges in a
wide variety of fields, including active matter [52, 53], air pollution [54],
astrophysics [55], biology [56, 57], chemistry [58, 59, 60, 61], cognitive
science [62], discrepancy modeling [63, 64, 65], epidemiology [66], fluid
dynamics [67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80],
manufacturing [81], heat transfer [82], numerical analysis [83], optics [84],
plasma physics [85, 86, 87], power systems [88, 89], structural dynamics [90,
91], and turbulence closure modeling [92, 93, 94, 95]. Furthermore, a large
number of innovations that extend the capabilities of the original SINDy
algorithm have been developed to discover partial differential equations
(PDEs) [50, 96], identify hybrid [97] and stochastic [98, 99, 100, 101, 102]
dynamical systems, incorporate control or exogenous inputs [103, 104, 105],
express system dynamics with rational functions [56, 106], enforce physical
constraints such as symmetries and conservation laws [67], global stability
[107], and dimensional consistency [108], add uncertainty quantification and
probabilistic forecasting capabilities [109, 110, 111, 112, 113], allow for
online model updating [114], leverage smart sampling strategies [115, 116],
propose model selection criteria [117, 118], and dramatically improve the
robustness of the method to noisy [119, 120, 121, 122, 123, 124], corrupt
[125, 126], and incomplete [127, 128, 129, 130] data. Many of these variants
have been implemented in the open source software package PySINDy [131, 132],
and some of them have been recently benchmarked [133] using the dysts
standardized database of chaotic systems introduced by [134]. Another
foundational contribution is that of Champion  _et al._ [135], who used
autoencoder neural networks along with SINDy to simultaneously discover latent
coordinates and dynamics from high-dimensional data. Subsequently, several
works have leveraged autoencoders in conjunction with SINDy for dimensionality
reduction and identifying dynamics, respectively, in novel methods for data-
driven low-order modeling for nonlinear PDEs [73, 136, 77, 137].

In the context of parameterized systems, the original SINDy article already
introduced a basic method that is able to successfully identify parameterized
dynamics, if provided with a large and noise-free data set [1]. Recently,
Meidani and Farimani [138] used integer programming to drastically improve the
noise robustness of the method, albeit for the specific scenario where the
model coefficients are known to be integers. The works of Schaeffer  _et al._
[139] and Rudy  _et al._ [140] used group sparsity to identify dynamical
systems and PDEs, respectively, from data sampled at multiple parameter
values. The group sparsity results in the identification of one model
structure with model coefficients that are allowed to change as the parameter
of the system varies, thus enabling model discovery in a scenario where the
parameter is unknown, slowly drifting, or abruptly changing [139, 140].
Consequently, it is important to point out that, this approach does not
produce a symbolic expression for the model dependence on the parameter. In
recent work, Nicolaou  _et al._ [141] developed a simple extension of the PDE-
FIND algorithm [50] to identify parameterized PDEs by carefully building the
library of function candidates to represent the parameter dependence. This
extension was shown to successfully produce models capable of extrapolation to
predict bifurcations [141]. However, the authors observe that the main
challenges to discover accurate models with this approach are that the method
requires sufficiently informative trajectory data and that sufficiently many
parameter values should be measured [141]. Interestingly, they also remark
that, according to their findings, trajectories with persistent dynamics
provide more information than transient trajectories. Also recently, Conti
_et al._ [137] developed a data-driven method to build nonlinear low-order
models of parameterized PDEs. They leveraged the autoencoder + SINDy framework
[135, 73, 136, 77] to identify parameterized dynamics on latent coordinates in
a low-dimensional embedding [137], which allowed for computationally cheap
numerical continuation of periodic solutions. Importantly, successful
identification of parameterized dynamical systems with any of the methods
described above requires a rich dataset containing measurements of multiple
transient trajectories with samples distributed across parameter space.
Reducing the amount of transient data required for parameterized model
discovery is the main subject of this work.

Figure 1: Schematic of the multi-objective SINDy framework for parameterized
systems. The method balances measurements on and off attractors to accurately
identify parameterized dynamics from sparse data. In the example, a Hopf
normal form is correctly identified from a single transient trajectory
complemented with a few measurements of stable fixed points and limit cycles
distributed over parameter space.

In an experimental setting, obtaining data on transient trajectories amounts
to acquiring time-resolved measurements, whereas data on attractors is
obtained from long-term response measurements. In addition, a stable fixed
point can be accurately measured using non-time-resolved techniques, and a
stable periodic orbit can be observed using stroboscopic measurements or phase
averaging time-resolved data. In a numerical setting, we generate transient
trajectories by time-marching through forward integration, whereas we may use
iterative solvers to find fixed points and periodic orbits. Moreover,
numerical continuation can be used to efficiently follow these solution
branches through parameter space [142]. In short, different experimental
techniques and numerical methods are typically used to obtain data on or off
attractors. As a very general statement, measurements on attractors can
typically be acquired at a lower computational cost (or with a simpler
experimental setup) and often have a higher fidelity, compared to transient
recordings of the same system. Furthermore, attractors and their bifurcations
form the backbone of the state space of nonlinear dynamical systems and,
because of this, we hypothesize that they provide more valuable data samples
for parameterized system identification.

In this work, we develop a SINDy extension to identify parameterized dynamical
systems using mostly data on attractors, and a minimum amount of data from
transient trajectories. Our proposed approach, schematically shown in figure
1, introduces a multi-objective sparse regression that balances the fit to
both types of measurements to drastically improve the numerical conditioning
of the regression.

The remainder of the paper is organized as follows. Theoretical background on
SINDy for parameterized systems and SINDy with constraints is covered in §II.
Our multi-objective SINDy framework is formulated in §III. Application of the
method to several numerical examples is presented and discussed in §V and our
conclusions are offered in §VI.

##  II Background

In this section we provide a brief review on the original SINDy formulation
for parameterized systems presented in Brunton  _et al._ [1] and the approach
to incorporate constraints introduced in Loiseau and Brunton [67].

###  II.1 SINDy for parameterized systems

The starting point is a parameterized nonlinear dynamical system with unknown
dynamics that we want to identify. SINDy assumes that the dynamics can be
expressed as a linear combination of a few nonlinear functions selected among
a much larger, and judiciously designed, set of candidate functions, as
follows

$\bm{\dot{x}}=\bm{f}(\bm{x},\mu)\approx\bm{\Xi}^{\mathrm{T}}\bm{\theta}(\bm{x},\mu),$
(1)

where the overdot denotes time-differentiation, $\bm{x}\in\mathbb{R}^{n}$ is
the state of the system, $\mu\in\mathbb{R}$ is a parameter,
$\bm{f}:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}^{n}$ corresponds
to the unknown dynamics,
$\bm{\theta}:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}^{d}$ contains
the set of $d$ candidate scalar functions, and $\bm{\Xi}\in\mathbb{R}^{d\times
n}$ are the yet to be identified model coefficients.

Given a set of $m$ measurements of the state and parameter
$\\{\bm{x}_{j},\mu_{j}\\}$ for $j=1,\dots,m$, typically acquired from several
transient trajectories generated for different parameter values, we may
assemble the data matrices

$\bm{\dot{X}}=\left[\bm{\dot{x}}_{1}\ \bm{\dot{x}}_{2}\ \cdots\
\bm{\dot{x}}_{m}\right]^{\mathrm{T}}\in\mathbb{R}^{m\times n},\quad\text{and
}\quad\bm{\Theta}=\left[\bm{\theta}(\bm{x}_{1},\mu_{1})\
\bm{\theta}(\bm{x}_{2},\mu_{2})\ \cdots\
\bm{\theta}(\bm{x}_{m},\mu_{m})\right]^{\mathrm{T}}\in\mathbb{R}^{m\times d},$

where the time derivatives $\bm{\dot{x}}_{j}=\dot{\bm{x}}(t_{j})$ may be
approximated from sequential state data, via finite differences for example,
if not directly measured. We may now formulate an optimization problem to
identify $\bm{\Xi}$, as follows

$\bm{\Xi}=\underset{\bm{\Xi^{\prime}}}{\mathrm{argmin}}\
\|\bm{\dot{X}}-\bm{\Theta}{\bm{\Xi^{\prime}}}\|_{F}^{2}+\lambda\|\bm{\Xi^{\prime}}\|_{0},$
(2)

where the first term in the cost function accounts for how well the model fits
the data, and the second term, weighted by the hyperparameter $\lambda$,
penalizes the amount of non-zero entries in the coefficient matrix, that is,
the amount of terms in the identified model. Finding an exact solution to (2)
requires performing a brute force search over all possible model structures,
which is intractable when the number of candidate functions is large. In
practice, sparse regression techniques, such as LASSO [51] or the sequentially
thresholded least squares (STLS) introduced in the original SINDy article [1],
bypass the combinatorial search by only approximating the optimal solution to
(2).

Note that this formulation slightly differs from that presented for
parameterized systems in the original article. In [1] the state is simply
augmented to contain the parameter as a new state variable governed by the
trivial dynamics $\dot{\mu}=0$, which is included as an additional equation
for the system. Although this may be convenient from a programming standpoint,
because it allows reusing the same code developed for non-parameterized
systems, it is not strictly necessary and the extra equation may be avoided
altogether.

It is important to point out once again that robust identification of a
parameterized dynamical system with SINDy typically requires a rich dataset
comprised of multiple transient trajectories generated for multiple parameter
values over the range of interest. In this work we develop a method to
significantly reduce the amount of transient data required for successful
identification of parameterized systems.

###  II.2 SINDy with constraints

The method introduced by Loiseau and Brunton [67] allows incorporating
physical priors in the identification process by solving a constrained sparse
regression. This is achieved by modifying the optimization problem (2) to

$\bm{\xi}=\underset{\bm{\xi^{\prime}}}{\mathrm{argmin}}\
\|\bm{\dot{\chi}}-\bm{\hat{\Theta}}{\bm{\xi^{\prime}}}\|_{2}^{2}+\lambda\|\bm{\xi^{\prime}}\|_{0},\quad\text{subject
to}\quad\bm{C}\bm{\xi^{\prime}}=\bm{d},$ (3)

where $\bm{\xi}$ and $\bm{\dot{\chi}}\in\mathbb{R}^{nd}$ are vectorized forms
of the $\bm{\Xi}$ and $\bm{\dot{X}}$ matrices, respectively,
$\bm{\hat{\Theta}}\in\mathbb{R}^{nm\times nd}$ is a block-diagonal matrix with
$n$ copies of $\bm{\Theta}$ as its blocks, and the matrix
$\bm{C}\in\mathbb{R}^{p\times nd}$ and vector $\bm{d}\in\mathbb{R}^{p}$
enforce $p$ linear equality constraints. The optimization problem (3) is then
approximately solved by iteratively solving a constrained least squares
problem, finding small coefficients that satisfy $|\xi_{i}|<\lambda$, and
augmenting $\bm{C}$ and $\bm{d}$ to incorporate the constraint
$\bm{e}_{i}^{\mathrm{T}}\bm{\xi}=\xi_{i}=0$ to promote sparsity.

The constrained least squares regression obtained for each iteration is a
convex optimization problem, has an analytical solution that can be explicitly
built for small-scale problems, and can be solved using readily available
software for larger-scale ones [67]. Crucially, this solution exists if and
only if the following conditions are satisfied:

  1. 1.

$\bm{C}$ has linearly independent rows, and

  2. 2.

$\begin{bmatrix}\bm{\hat{\Theta}}\\\ \bm{C}\end{bmatrix}$ has linearly
independent columns.

The first condition requires that $\bm{C}$ is wide ($p\leq nd$) meaning that
there are fewer constraints than model coefficients. The second condition
depends on both, $\bm{\hat{\Theta}}$ and $\bm{C}$, and it can be satisfied
even if the function library has linearly dependent columns.

In addition to enforcing sparsity, constraints allow the incorporation of
problem-specific knowledge of the system, such as fixed points, symmetries,
and conservation laws [67]. The resulting models not only respect the known
physics, but also generalize better and require less data for training.
Moreover, in the context of identifying parameterized systems from as few
measurements as possible across parameter space, SINDy with constraints seems
like an attractive alternative to enforce parameter-dependent priors, such as
fixed points and/or limit cycles. However, as we discuss in the next section,
there are critical pitfalls to this approach.

##  III Proposed method

In this section, we develop our SINDy extension to identify parameterized
dynamical systems from limited data. We begin with our motivation to
distinguish between data obtained on and off attractors, continue with the
description of a first method that has limited use in practice, and then pivot
to our proposed multi-objective sparse regression.

###  III.1 Data on and off attractors

Acquisition of time-resolved measurements from transient trajectories and
long-term response measurements from attractors often require different
experimental setups. Even for the case of numerical data, we may solve for a
stable steady or periodic state without having to perform a transient
simulation. Furthermore, numerical continuation allows us to efficiently
follow attractors in parameter space. This translates into a significant
difference in cost, either economic or computational, between the acquisition
of data on and off attractors. Moreover, attractors and their bifurcations
provide a skeleton for the state space dynamics of parameterized systems. For
these reasons, we hypothesize that an accurate representation of attractors is
critical to successfully identify the correct model structure and, more so, it
may even be sufficient to capture the parameter-dependence of the dynamics.
This motivates the development of a method to identify parameterized dynamical
systems relying as much as possible on data on attractors, and using as less
data as possible from transient trajectories. Our approach to achieve this is
to enforce the on-attractor dynamics as constraints.

We begin by distinguishing between data from transient trajectories
$\\{\bm{x}^{\mathrm{tr}}_{j},\mu^{\mathrm{tr}}_{j}\\}$, for $j=1,\dots,m$, and
data on attractors $\\{\bm{x}^{\mathrm{att}}_{j},\mu^{\mathrm{att}}_{j}\\}$,
for $j=1,\dots,r$, leading to the definition of the following data matrices

$\displaystyle\bm{\dot{X}}_{\mathrm{tr}}=\left[\bm{\dot{x}}^{\mathrm{tr}}_{1}\
\bm{\dot{x}}^{\mathrm{tr}}_{2}\ \cdots\
\bm{\dot{x}}^{\mathrm{tr}}_{m}\right]^{\mathrm{T}},\quad$
$\displaystyle\bm{\Theta}_{\mathrm{tr}}=\left[\bm{\theta}(\bm{x}^{\mathrm{tr}}_{1},\mu^{\mathrm{tr}}_{1})\
\bm{\theta}(\bm{x}^{\mathrm{tr}}_{2},\mu^{\mathrm{tr}}_{2})\ \cdots\
\bm{\theta}(\bm{x}^{\mathrm{tr}}_{m},\mu^{\mathrm{tr}}_{m})\right]^{\mathrm{T}},$
$\displaystyle\bm{\dot{X}}_{\mathrm{att}}=\left[\bm{\dot{x}}^{\mathrm{att}}_{1}\
\bm{\dot{x}}^{\mathrm{att}}_{2}\ \cdots\
\bm{\dot{x}}^{\mathrm{att}}_{r}\right]^{\mathrm{T}},\quad\text{and}$
$\displaystyle\bm{\Theta}_{\mathrm{att}}=\left[\bm{\theta}(\bm{x}^{\mathrm{att}}_{1},\mu^{\mathrm{att}}_{1})\
\bm{\theta}(\bm{x}^{\mathrm{att}}_{2},\mu^{\mathrm{att}}_{2})\ \cdots\
\bm{\theta}(\bm{x}^{\mathrm{att}}_{r},\mu^{\mathrm{att}}_{r})\right]^{\mathrm{T}}.$

Furthermore, if the attractor is a fixed point, then the time derivative is
$\bm{\dot{x}}^{\mathrm{att}}_{j}=\bm{0}$, and, if it refers to a limit cycle,
then it can be computed analytically from a Fourier series representation.

###  III.2 The limitation of hard constraints

As a first effort to formulate a method that leverages data on attractors
across parameter space, we attempt using the constrained sparse identification
developed in [67] and discussed in §II. We start from the optimization problem
(3) using the data from the transient trajectories. In addition to the
constraints $\bm{C}\bm{\xi}=\bm{d}$, used to enforce the sparsity and other
knowledge about the system, we include additional linear equality constraints
$\bm{\hat{\Theta}}_{\mathrm{att}}\bm{\xi}=\bm{\dot{\chi}}_{\mathrm{att}}$ to
impose the on-attractor dynamics, as follows

$\bm{\xi}=\underset{\bm{\xi^{\prime}}}{\mathrm{argmin}}\
\|\bm{\dot{\chi}}_{\mathrm{tr}}-\bm{\hat{\Theta}}_{\mathrm{tr}}{\bm{\xi^{\prime}}}\|_{2}^{2}+\lambda\|\bm{\xi^{\prime}}\|_{0},\quad\text{subject
to}\quad\bm{\hat{\Theta}}_{\mathrm{att}}\bm{\xi^{\prime}}=\bm{\dot{\chi}}_{\mathrm{att}}\
\text{ and }\ \bm{C}\bm{\xi^{\prime}}=\bm{d},$ (4)

where $\bm{\hat{\Theta}}_{\mathrm{tr}}\in\mathbb{R}^{nm\times nd}$ and
$\bm{\hat{\Theta}}_{\mathrm{tr}}\in\mathbb{R}^{nr\times nd}$ are the block-
diagonal matrices with $n$ copies of the corresponding data matrices for the
function libraries, and $\bm{\dot{\chi}}_{\mathrm{tr}}\in\mathbb{R}^{nm}$ and
$\bm{\dot{\chi}}_{\mathrm{att}}\in\mathbb{R}^{nr}$ are the vectorized data
matrices with the time derivatives. Both sets of equality constraints can be
combined to obtain

$\bm{\hat{\Theta}}_{\mathrm{att}}\bm{\xi}=\bm{\dot{\chi}}_{\mathrm{att}}\
\text{ and }\
\bm{C}\bm{\xi^{\prime}}=\bm{d}\qquad\Leftrightarrow\qquad\bm{\tilde{C}}\bm{\xi}=\bm{\tilde{d}},\quad\bm{\tilde{C}}=\begin{bmatrix}\bm{\hat{\Theta}}_{\mathrm{att}}\\\
\bm{C}\end{bmatrix},\quad\bm{\tilde{d}}=\begin{bmatrix}\bm{\dot{\chi}}_{\mathrm{att}}\\\
\bm{d}\end{bmatrix}.$

Importantly, in order for a solution to the optimization problem to exist, we
require that the matrix $\bm{\tilde{C}}\in\mathbb{R}^{(p+nr)\times nd}$ is
wide. Consequently, the dimensions of $\bm{\tilde{C}}$ must satisfy
$(p+nr)\leq nd$, or equivalently,

$nr\leq nd-p,$ (5)

where $nr$ is the total number of constraints imposed to enforce the on-
attractor dynamics, $nd$ is the total number of model coefficients being fit,
and $p$ is the number of sparsity constraints, making $nd-p$ the amount of
non-zero model coefficients. Therefore, the number of constraints that one can
impose is, at most, the amount of active coefficients in the underlying
dynamics. However, the assumption that the dynamics are governed by equations
with only a few terms is at the core of the SINDy approach. Moreover, when
attempting to enforce the on-attractor dynamics across a parameter range, even
if it is a single fixed point measured at several parameter values, $nr$ grows
with the number of samples $r$.

As a consequence, condition (5) will rarely be satisfied in practice, making
this approach extremely limited in its application to enforce the
parameterized on-attractor dynamics. Therefore, we do not recommend using the
method presented in this subsection for this purpose. Instead, we propose the
approach presented next.

###  III.3 Multi-objective sparse regression

In view of the limited applicability of the method described above, we
consider a different approach. Rather than incorporating the on-attractor
dynamics as hard constraints, we include them in our cost function as soft
constraints that merely penalize, as opposed to strictly enforcing, the model
error on the attractor data. This leads to the formulation of a multi-
objective sparse regression problem, as follows

$\bm{\Xi}=\underset{\bm{\Xi^{\prime}}}{\mathrm{argmin}}\
\underbrace{\|\bm{\dot{X}}_{\mathrm{tr}}-\bm{\Theta}_{\mathrm{tr}}{\bm{\Xi^{\prime}}}\|_{F}^{2}}_{\textstyle{J_{\mathrm{tr}}(\bm{\Xi^{\prime}})}}+\alpha\underbrace{\|\bm{\dot{X}}_{\mathrm{att}}-\bm{\Theta}_{\mathrm{att}}{\bm{\Xi^{\prime}}}\|_{F}^{2}}_{\textstyle{J_{\mathrm{att}}(\bm{\Xi^{\prime}})}}+\lambda\|\bm{\Xi^{\prime}}\|_{0},$
(6)

where the hyperparameter $\alpha$ weights how much we emphasize on minimizing
the fit error on the attractors $J_{\mathrm{att}}(\bm{\Xi})$ compared to the
error on the transient trajectories $J_{\mathrm{tr}}(\bm{\Xi})$. Note that,
for a very small value of $\alpha$, the optimization will disregard the data
on the attractors, whereas, for a very large $\alpha$, it will deemmphasize
the transients. In fact, in the limit as $\alpha\rightarrow\infty$ we recover
the constrained sparse regression (4), which is ill posed if (5) is not
satisfied. Therefore, as we show in our numerical examples in §V, there is
typically a sweet spot for $\alpha$ that balances the fit to attractors and
transients to find an accurate model across parameter space. As with any other
hyperparameter, $\alpha$ may be selected via cross-validation or more
sophisticated model selection criteria [117, 118].

Furthermore, the cost function in (6) may be rewritten to yield

$\bm{\Xi}=\underset{\bm{\Xi^{\prime}}}{\mathrm{argmin}}\
\left\|\vphantom{\begin{bmatrix}.\\\
.\end{bmatrix}}\right.\underbrace{\begin{bmatrix}\bm{\dot{X}}_{\mathrm{tr}}\\\
\sqrt{\alpha}\bm{\dot{X}}_{\mathrm{att}}\end{bmatrix}}_{\textstyle{\bm{\dot{X}}}}-\underbrace{\begin{bmatrix}\bm{\Theta}_{\mathrm{tr}}\\\
\sqrt{\alpha}\bm{\Theta}_{\mathrm{att}}\end{bmatrix}}_{\textstyle{\bm{\Theta}}}{\bm{\Xi^{\prime}}}\left.\vphantom{\begin{bmatrix}.\\\
.\end{bmatrix}}\right\|_{F}^{2}+\lambda\|\bm{\Xi^{\prime}}\|_{0},$ (7)

which has an equivalent structure to problem (2) and, for the case $\alpha=1$,
we recover the original SINDy formulation [1] using all the measurements
stacked together. Thus, for a given $\alpha$, we can approximate the solution
to (7) via STLS. In fact, we me reuse any available SINDy code after forming
the row-weighted data matrices $\bm{\dot{X}}$ and $\bm{\Theta}$ in (7). The
pseudocode for our implementation is shown in algorithm 1, where we use
Matlab’s array slicing notation.

Algorithm 1 Multi-objective SINDy for parameterized systems

Inputs: data from transients $\bm{\Theta}_{\mathrm{tr}}$ and
$\bm{\dot{X}}_{\mathrm{tr}}$

data on attractors $\bm{\Theta}_{\mathrm{att}}$ and
$\bm{\dot{X}}_{\mathrm{att}}$

sparsification and trade-off hyperparameters $\lambda$ and $\alpha$

Outputs: SINDy model coefficients $\bm{\Xi}$

1:$\bm{\dot{X}}\leftarrow[\bm{\dot{X}}_{\mathrm{tr}};\sqrt{\alpha}\bm{\dot{X}}_{\mathrm{att}}]$
$\triangleright$ Weight and vertically concatenate the time derivatives.

2:$\bm{\Theta}\leftarrow[\bm{\Theta}_{\mathrm{tr}};\sqrt{\alpha}\bm{\Theta}_{\mathrm{att}}]$
$\triangleright$ Weight and vertically concatenate the function libraries.

3:function STLS($\bm{\dot{X}}$,$\bm{\Theta}$,$\lambda$)

4: $\bm{v}\leftarrow\max(\mathrm{abs}(\mathrm{columns}(\bm{\Theta})))$
$\triangleright$ Find the largest entry for each candidate function…

5: $\bm{\Theta}\leftarrow\bm{\Theta}/\mathrm{diag}(\bm{v})$ $\triangleright$
…and normalize the columns of the function library.

6: $\bm{\Xi}\leftarrow\bm{\Theta}\backslash\bm{\dot{X}}$ $\triangleright$
Initialize the coefficients with the least squares solution.

7: for $k=1,\dots,5$ do $\triangleright$ During each iteration…

8: for $i=1:n$ do $\triangleright$ …and for each state variable…

9:
$\bm{j}_{\mathrm{small}}\leftarrow\mathrm{abs}(\bm{\Xi}(:,i))<\lambda\max(\mathrm{abs}(\bm{\Xi}(:,i)))$
$\triangleright$ …find the entries with a small contribution to the dynamics…

10: $\bm{\Xi}(\bm{j}_{\mathrm{small}},i)\leftarrow\bm{0}$ $\triangleright$
…and threshold.

11: $\bm{j}_{\mathrm{big}}\leftarrow\ \sim\bm{j}_{\mathrm{small}}$
$\triangleright$ Then find the non-zero entries…

12:
$\bm{\Xi}(\bm{j}_{\mathrm{big}},i)\leftarrow\bm{\Theta}(:,\bm{j}_{\mathrm{big}})\backslash\bm{\dot{X}}(:,i)$
$\triangleright$ …and regress onto those terms.

13: end for

14: end for

15: $\bm{\Xi}\leftarrow\mathrm{diag}(\bm{v})\backslash\bm{\Xi}$
$\triangleright$ Rescale the rows of the coefficent matrix.

16:end function

There are two subtle details in the implementation of the STLS procedure that
we use in algorithm 1. First, because we have row-weighted data matrices,
normalization of the columns of $\bm{\Theta}$ becomes necessary to avoid
having to scale the sparsity knob $\lambda$ with $\alpha$. Second, we
introduce a variable-dependent thresholding criterion, where we retain the
terms that most contribute to the dynamics of each state variable, rather than
using one global thresholding criterion. We find that this becomes relevant
when there is a significant variation in the size of the coefficients between
the equations for different state variables.

A schematic of the application of the method to identify a Hopf normal form is
shown in figure 1. Our multi-objective SINDy framework successfully identifies
the parameterized dynamics from data of a single transient trajectory along
with stable fixed point and limit cycle data acquired for several parameter
values. Due to its simplicity, the proposed approach may be trivially coupled
with other powerful SINDy extensions, for example to identify PDEs [50], add
noise robustness [96], and quantify model uncertainty [112].

##  IV Numerical examples and dataset

To demonstrate the application of multi-objective SINDy, we generate a dataset
comprised of numerical measurements of four parameterized dynamical systems: a
saddle-node hysteresis loop, the normal form of a supercritical Hopf
bifurcation, a system of coupled Stuart–Landau oscillators, and the Lorenz
equations. These systems are selected to feature all kinds of attractors;
fixed points, limit cycles, quasiperiodic orbits, and strange attractors. For
every system, we generate data by choosing an initial condition and
numerically integrating long trajectories for a range of parameter values.
Numerical integration is carried out using an explicit $4-5$th order
Runge–Kutta method implemented in the SciPy library in Python. Data points on
and off attractors are registered, and a single transient trajectory,
highlighted in figure 2, is flagged to be used for parameterized model
discovery in §V. Data on transient trajectories are acquired at a fast
sampling rate to appropriately resolve the time evolution. Time derivatives of
the state variables are approximated with a $4$th order central finite
difference scheme. Data on limit cycles, quasiperiodic and chaotic orbits are
obtained using a burst sampling strategy [116], where short bursts of
measurements are acquired at a fast sampling rate (the same used for
transients) and separated by a much slower sampling rate to provide a better
coverage of the attractor. Only the central data point within each burst is
retained and their neighbors are used to approximate the corresponding time
derivative and then discarded. Details on the selected dynamical systems and
their data are presented below.

Figure 2: Numerical dataset of parameterized dynamical systems used to develop
and compare methods for parameterized model discovery. For each system, off-
and on-attractor measurements are acquired over a range of parameter values.
The highlighted single transient trajectories are used for parameterized model
discovery in §V. This dataset is available at github.com/ben-herrmann.

###  IV.1 Saddle-node hysteresis loop

The dynamics for the saddle-node hysteresis loop are given by

$\dot{x}=\mu+x-x^{3}.$ (8)

Qualitatively similar models have been used to describe the dynamics of
thermal runaway in batteries and insect outbreaks in forests [143]. This
system exhibits a pair of saddle-node bifurcations at $\mu=\pm 1/\sqrt{3}$,
leading to two different branches of stable equilibria that coexist in the
range $-1/\sqrt{3}<\mu<1/\sqrt{3}$.

Data for this system is acquired for parameter values starting at $\mu=-6$ and
reaching $\mu=-0.182$ in steps of $0.65$, and also starting at $\mu=0.182$ and
reaching $\mu=6$ in steps of $0.65$, for a total of $20$ sampled parameter
values. For each $\mu$, transient trajectories are integrated starting from
the initial condition $x_{0}=0$ and $1000$ data points are registered using a
sampling time of $0.01$ time units. A single data point is registered on one
of the stable equilibria branches for every parameter value. The negative
branch is sampled for negative $\mu$ values, and the positive one for positive
$\mu$ values. The single transient selected to test the proposed method
corresponds to $\mu=-0.182$.

###  IV.2 Hopf normal form

The governing equations for a supercritical Hopf normal form are

$\displaystyle\dot{x}$ $\displaystyle=\mu x-y-x(x^{2}+y^{2}),$ (9a)
$\displaystyle\dot{y}$ $\displaystyle=x+\mu y-y(x^{2}+y^{2}).$ (9b)

This is a well-known model for ubiquitous self-sustained oscillations
appearing in fluid dynamics, chemistry and biology. The bifurcation occurs at
$\mu=0$, where the fixed point $(0,0)^{\mathrm{T}}$ looses its stability and
trajectories end up in a stable limit cycle given by
$\bm{x}(t)=\sqrt{\mu}(\sin(t+\phi),\cos(t+\phi))^{\mathrm{T}}$ with a phase
$\phi$ that depends on the initial conditions.

Data for this system is acquired for parameter values starting at $\mu=-3$ and
increasing up to $\mu=2$ in steps of $0.5$, and also at $\mu=0.3,\ 0.4,$ and
$0.5$, for a total of $14$ sampled values. The value $\mu=0$ lands exactly at
the bifurcation point so it is replaced by $\mu=0.01$. For each $\mu$,
transient trajectories are integrated starting from the initial condition
$\bm{x}_{0}=(2,0)^{\mathrm{T}}$ and $2000$ data points are registered using a
sampling time of $0.01$ time units. A single data point is registered on the
stable equilibrium for each $\mu<0$, and, after the bifurcation, $5$ data
points are acquired on the limit cycle sampled every $1.1$ time units for each
$\mu>0$. The single transient selected to produce the results in §V
corresponds to $\mu=0.01$.

###  IV.3 Coupled Stuart–Landau oscillators

We consider two coupled Stuart–Landau oscillators governed by

$\displaystyle\dot{x}_{1}$
$\displaystyle=-\omega_{1}y_{1}+\left(R_{1}^{2}-(x_{1}^{2}+y_{1}^{2})\right)x_{1}+\mu(x_{1}x_{2}+y_{1}y_{2}),$
(10a) $\displaystyle\dot{y}_{1}$
$\displaystyle=\omega_{1}x_{1}+\left(R_{1}^{2}-(x_{1}^{2}+y_{1}^{2})\right)y_{1}+\mu(x_{1}y_{2}-y_{1}x_{2}),$
(10b) $\displaystyle\dot{x}_{2}$
$\displaystyle=-\omega_{2}y_{2}+\left(R_{2}^{2}-(x_{2}^{2}+y_{2}^{2})\right)x_{2},$
(10c) $\displaystyle\dot{y}_{2}$
$\displaystyle=\omega_{2}x_{2}+\left(R_{2}^{2}-(x_{2}^{2}+y_{2}^{2})\right)y_{2}.$
(10d)

A similar set of equations were recently proposed as a universal model for
periodically forced turbulent oscillator flows [144]. Here, the unforced
frequencies and amplitudes of each oscillator are fixed to $R_{1}=1,\
R_{2}=0.2,\ \omega_{1}=1/\pi,$ and $\omega_{2}=1$, respectively. The only free
parameter left is $\mu$ that represents the strength of the coupling between
the oscillators. Note that the coupling is only in one direction, with the
second oscillator driving the first one. When the second oscillator reaches
its periodic orbit, this driving acts as periodic forcing on the first
oscillator. Because the frequencies $\omega_{1}$ and $\omega_{2}$ are
incommensurate, if the forcing is weak, the system evolves towards a
quasiperiodic attractor. At $\mu=|\omega_{2}-2\omega_{1}|/(2R_{2})\approx
0.909$ the system undergoes a bifurcation and, for stronger forcing, a
phenomenon known as $2:1$ synchronization emerges [144]. What happens is that
the first oscillator phase locks to half the frequency of the second one,
leading to a periodic solution corresponding to a stable limit cycle that
appears as the underlying attractor.

Data for this system is acquired for $10$ parameter values ranging from
$\mu=0$ to $\mu=1.817$ in steps of $0.202$. For each $\mu$, transient
trajectories are integrated starting from the initial condition
$\bm{x}_{0}=(10^{-4},0,10^{-2},0)^{\mathrm{T}}$ and $5020$ data points are
registered using a sampling time of $0.01$ time units. After integrating for
$2000$ time units, a trajectory is considered to be on the underlying
attractor, and $12$ data points sampled every $1$ time units are registered
for each $\mu$. The single transient selected to produce the results in §V
corresponds to $\mu=1.817$.

###  IV.4 Lorenz system

We consider the famous Lorenz equations

$\displaystyle\dot{x}$ $\displaystyle=\sigma(y-x),$ (11a)
$\displaystyle\dot{y}$ $\displaystyle=x(\mu-z)-y,$ (11b)
$\displaystyle\dot{z}$ $\displaystyle=xy-\beta z,$ (11c)

that represent a simple model of heat convection in the atmosphere [145]. We
take $\sigma=10$, $\beta=8/3$, and $\mu$ takes the place of the parameter
classically denoted as $\rho$. We focus in the range $16\leq\mu\leq 28$, where
the system has $2$ stable equilibria at
$(\pm\sqrt{\beta(\mu-1)},\pm\sqrt{\beta(\mu-1)},\mu-1)^{\mathrm{T}}$ that
become unstable after undergoing a subcritical Hopf bifurcation when
increasing the parameter over $\mu\approx 24.74$, leading to the emergence of
a chaotic or strange attractor.

Data for this system is acquired for $\mu=16$ up to $\mu=28$ in increments of
$1$, for a total of $12$ sampled parameter values. For each $\mu$, transient
trajectories are integrated starting from the initial condition
$\bm{x}_{0}=(30,-30,50)^{\mathrm{T}}$ and $500$ data points are registered
using a sampling time of $0.001$ time units. For $\mu<24.74$, one data point
is registered for only on of the stable fixed points. For values of the
parameter after the bifurcation, we integrate long trajectories and, after
$70$ time units, we consider them to be on the chaotic attractor and register
$25$ data points sampled every $0.2$ time units. The transient trajectory
corresponding to $\mu=26$ is selected to test the multi-objective SINDy method
in §V.

##  V Results and discussion

Figure 3: Multi-objective SINDy applied to a number of examples. Results are
obtained using data from a single transient trajectory and from attractors
measured over a parameter range. Performance at parameterized model discovery
is assessed in terms of model coefficient error and the correctness of the
identified model structure. Multi-objective SINDy models perform significantly
better than standard SINDy models trained on the same data, and almost as well
as standard SINDy models trained on a much larger dataset. This is explained
because the row-weighting provided by the trade-off hyperparameter $\alpha$
drastically improves the numerical conditioning of the function library data
matrices.

In this section, we demonstrate the application of multi-objective SINDy using
the data for the systems described in the previous section. First, we apply
our proposed method to identify the parameterized dynamics of these systems
using a single transient trajectory and data on the corresponding attractors,
shown in figure 3. We then investigate the noise robustness and data
requirements of the method, shown in figure 4.

For each system, time derivatives are approximated via finite differences
using a $4$th order central scheme. As mentioned in the previous section, time
derivatives for data points on the attractors rely on the burst sampling
strategy [116]. For the function libraries, we use monomials of up to degrees
$5,\ 4,\ 3$, and $4$ in the state variables and parameters for the saddle-node
hysteresis loop, Hopf normal form, coupled Stuart-Landau oscillators, and
Lorenz system, respectively.

###  V.1 Parameterized system identification from a single transient
trajectory

The model identification performance of our method is compared to that of the
standard SINDy algorithm from [1] using a combination of data on the
attractors and from a single transient trajectory, as shown in figure 3.
Furthermore, to provide a more complete picture, we also present the
performance of the standard SINDy algorithm using only the data on the
attractors, only the data from the single transient trajectory, and the full
dataset, shown in figure 2, containing multiple transient trajectories. Only
noise-free numerical data is considered in this section. Models are identified
for each system and dataset using various values for the hyperparameters
$\lambda$ and $\alpha$. For the the sparsification knob $\lambda$ we do a
sweep over a grid of $100$ log-spaced values ranging from $10^{-4}$ to $1$.
When using multi-objective SINDy, the trade-off hyperparameter $\alpha$ is
evaluated over a grid of $11$ log-spaced values between $10^{-20}$ and
$10^{20}$. The hyperparameters that produce the best performing models for
each case are shown in table 1.

System |  $\lambda$ (SINDy) |  $\lambda,\ \alpha$ (Multi-obj.)
---|---|---
Single transient | Attractors | Single trans.$+$attr. | Full dataset | Single trans.$+$attr.
Saddle-node hysteresis loop | $1.000$ | $0.0001$ | $0.0105$ | $0.0095$ | $0.0045,\ 10^{20}$
Hopf normal form | $1.000$ | $0.3944$ | $0.2719$ | $0.0003$ | $0.0060,\ 10^{4}$
Coupled Stuart–Landau oscillators | $0.0001$ | $0.0152$ | $0.0034$ | $0.0001$ | $0.0087,\ 10^{4}$
Lorenz system | $0.0072$ | $0.0001$ | $0.0001$ | $0.0001$ | $0.0001,\ 10^{20}$
Table 1: Hyperparameters used to compare the performance of SINDy and multi-
objective SINDy. For each system and training data, these values produce the
lowest model coefficient error over a hyperparameter sweep.

The performance of the identified models is assessed in terms of model
coefficient error defined as

$\text{Error}=\frac{\|\bm{\Xi}-\bm{\Xi}_{\mathrm{true}}\|_{F}}{\|\bm{\Xi}_{\mathrm{true}}\|_{F}},$
(12)

where $\bm{\Xi}$ are the identified model coefficients and
$\bm{\Xi}_{\mathrm{true}}$ are the true coefficients from the underlying
system. In each case, we assess whether or not the correct model structure was
identified based on the sparsity pattern of $\bm{\Xi}$. Results shown in
figure 3 correspond to the best performing models in each case using the
hyperparameters from table 1.

Notice that, in general, SINDy is not able to successfully identify a
parameterized dynamical system neither from a single transient trajectory nor
from data on attractors only. This is expected because, in the first case we
have not provided any information about the parameter dependence of the
system. In the second case, identification fails because SINDy finds the model
with the fewest number of terms to fit the data on the attractors, which is
often a simpler model than the one that also explains the off-attractor
dynamics. To provide some intuition, the simplest model that explains a fixed
point measurement is given by the trivial dynamics $\bm{\dot{x}}=\bm{0}$, and
the simplest model that explains observations from a limit cycle is given by a
linear model with a center equilibrium.

The exception to this is the Lorenz system that has a chaotic attractor that
provides enough information for successful for parameterized system
identification. However, for a system where the chaotic attractor lies on a
subspace of the full state space, measurements of transient trajectories would
still be required to identify the off-attractor dynamics. Such an example of a
chaotic system where data on the chaotic attractor is not enough for model
discovery can be easily produced. For example, this would be the scenario for
the Lorenz system if we added a fourth state variable $w$ with dynamics
$\dot{w}=-w$. On the attractor we would only measure $w=0$ and therefore would
not be able to identify the underlying dynamics.

Moreover, for the presented systems, standard SINDy is also not able to
correctly identify the underlying dynamics when using the combined data from a
single transient and attractors, and requires the full dataset containing
multiple transient trajectories, as shown in figure 3. We find that just by
weighting the data matrices, multi-objective SINDy correctly identifies the
model structure using data from just a single transient trajectory and on the
attractors. Furthermore, the performance in terms of model coefficient error
is almost as good as that obtained by the standard SINDy algorithm using the
full dataset.

Again, the Lorenz system is the exception to this general behavior.
Interestingly, for this case, data on the attractors is sufficient for
successful identification and, more so, including the transient trajectories
has a detrimental effect on model performance. Nevertheless, with multi-
objective SINDy we achieve a similar performance because the weighting by
$\alpha$ diminishes the relative contribution of the transient data to the
regression results.

Figure 4: Robustness to noise and data length requirements for parameterized
model discovery using data on attractors combined with sparse noisy data from
a single transient trajectory over a range of parameters. Performance metrics
are averaged over $500$ random realizations of measurement noise and data
removal. (a) Saddle-node hysteresis loop. (b) Hopf normal form. (c) Coupled
Stuart–Landau oscillators. (d) Lorenz system. Comparison between standard
SINDy and multi-objective SINDy shows that the proposed method is able to
identify more accurate models using less and more noisy data.

To further understand why the proposed method performs so well in comparison
the the standard approach, we compute the condition numbers of the function
library data matrices $\kappa(\bm{\Theta})$ obtained using different parts of
the dataset. In the last column of figure 3, we compare the condition numbers
of the data matrices formed using the full dataset
$\bm{\Theta}_{\mathrm{full}}$ (in black), data from a single transient
trajectory $\bm{\Theta}_{\mathrm{tr}}$ (in magenta), only on the attractors
$\bm{\Theta}_{\mathrm{att}}$ (in blue), and the combination of the latter two
weighted by the trade-off hyperparameter
$\bm{\Theta}_{\alpha}=\left[\bm{\Theta}_{\mathrm{tr}};\sqrt{\alpha}\bm{\Theta}_{\mathrm{att}}\right]$
(in purple) as a function of $\alpha$. The last is the data matrix involved in
the multi-objective sparse regression problem from equation (7) solved by the
proposed method, and corresponds to the data matrix used in standard SINDy
when $\alpha=1$. The plots for the Hopf normal form and coupled Stuart–Landau
oscillators in figure 3 clearly show that in the limit as
$\alpha\rightarrow-\infty$ the condition number
$\kappa(\bm{\Theta}_{\alpha})\rightarrow\kappa(\bm{\Theta}_{\mathrm{tr}})$, in
the limit as $\alpha\rightarrow\infty$ then
$\kappa(\bm{\Theta}_{\alpha})\rightarrow\kappa(\bm{\Theta}_{\mathrm{att}})$,
and that there is a minimum in between that gets much closer to
$\kappa(\bm{\Theta}_{\mathrm{full}})$ than either extrema. Importantly, this
minimum in the condition number does not necessarily occur at $\alpha=1$, and
its position will depend on the relative amount of data points from each type
(on- or off-attractor dynamics) and the data values themselves. As a
consequence, the multi-objective approach significantly improves the
conditioning of the regression problem by introducing a value of $\alpha\neq
1$ to better balance the information provided by the different types of
measurements. For the Lorenz system we find that the library evaluated on data
on the attractors has a condition number that is almost as low as that of the
library evaluated in the full dataset and is significantly better conditioned
than the library evaluated on a single transient. Moreover, for values of
$\alpha<1$, increasing the relative contribution of the transient
trajectories, we find that the conditioning drastically worsens, explaining
the detrimental effect of including the transient measurements on model
performance. Results shown here use the values of $\alpha$, shown in table 1,
that produce the best performing model in each case with the aim of
highlighting the potential of multi-objective SINDy. However, in a real system
identification scenario, we propose selecting the value of $\alpha$ via cross-
validation or, as with any other hyperparameter, it can also be determined
using other model selection criteria [117, 118].

###  V.2 Robustness to noise and data requirements

Improving the conditioning of the regression should provide robustness to the
identification procedure when dealing with sparse and noisy measurements.
Following the methodology of Fasel  _et al._ [112], we demonstrate this
robustness by adding random measurement noise to our data and assessing the
average performance over multiple realizations for a range of noise
intensities and data lengths. Specifically, for each of the systems presented
in §IV, we use the data from a single transient trajectory and on the
attractors for various values of the parameter. We add Gaussian white noise
with zero mean and variance $\sigma^{2}$ to the transient only for a range of
noise levels, defined as $\sigma$ normalized by the root mean square of the
transient trajectory data. In addition, we randomly discard measurements of
this same transient trajectory and retain a certain percentage of the original
number of data points. The model coefficient error and success rate ($1$ if
the correct model structure is identified, $0$ if not) are averaged over $500$
random realizations of noise and of discarded data points for each noise level
and percentage of data removal, as shown in figure 4. For each system, the
hyperparameters are fixed for all noise levels and data removal percentages to
the values shown in table 2. These values are lightly tuned to make the trends
in figure 4 clear enough to support our main conclusions. However, we remark
that these values have not been optimized in any way and that the performance
of multi-objective SINDy could probably be improved even more. Furthermore, in
a realistic scenario, the hyperparameters would be selected by cross-
validating against data with the same noise level and sparsity as the training
data, leading to hyperparameters that are specifically tuned for the
characteristics of the available measurements.

System |  $\lambda$ (SINDy) |  $\lambda,\ \alpha$ (Multi-objective SINDy)
---|---|---
Saddle-node hysteresis loop | $0.15$ | $0.15,\ 10^{14}$
Hopf normal form | $0.06$ | $0.06,\ 10^{1}4$
Coupled Stuart–Landau oscillators | $0.0087$ | $0.0087,\ 10^{1}0$
Lorenz system | $0.01$ | $0.01,\ 10^{14}$
Table 2: Hyperparameters used to compare the robustness to noise and data
requirements of SINDy and multi-objective SINDy. Models are trained using data
on attractors combined with sparse noisy data from a single transient
trajectory

As mentioned in the previous subsection, using a single transient trajectory
and the attractors, SINDy is not able to identify the parameterized dynamics
for the saddle-node hysteresis loop, Hopf normal form or coupled Stuart–Landau
oscillators, even when the data is clean and full length, that is with a noise
level of $0\%$ and $0\%$ data removal. With multi-objective SINDy, not only
can we identify these systems from clean full-length data, but the method
significantly moves the front of successful identification towards the top-
right corner of high intensity noise and sparse measurements, as shown in
figure 4. Furthermore, although SINDy is able to identify the parameterized
Lorenz dynamics without noise, multi-objective SINDy is more robust to noise.
As discussed in the previous section for the Lorenz example, including the
data on the transient trajectory, as opposed to just using data on the
attractors, is actually detrimental for system identification performance.
Notably, this effect is not seen for multi-objective SINDy, that deemphasizes
the data from the transient trajectory via the row-weighting of the data
matrices to achieve better performance.

##  VI Conclusions

Measurements of the state of a dynamical system evolving on an attractor
typically have a lower acquisition cost and are more informative for
parameterized system identification than measurements from transient
trajectories. In this work, we developed a method to identify parameterized
dynamical systems relying as much as possible on data on attractors, and using
as less data as possible from transient trajectories. The approach is an
extension of the SINDy method [1] that introduces a multi-objective cost
function to balance the fit to data on and off attractors, leading to a
significant improvement in the numerical conditioning of the underlying
regression problem. We show that parameterized nonlinear dynamics can be
learned from a single transient trajectory and data on attractors acquired
over a range of parameter values. The new multi-objective SINDy is more robust
to noisy measurements and requires less data for successful model discovery
than the standard SINDy algorithm.

Moreover, the developed method is very simple to implement and can be easily
coupled with other available SINDy extensions. Furthermore, the key idea
behind multi-objective SINDy can be adopted for problems beyond parameterized
dynamical systems, in other contexts where data from multiple sources can be
combined for model discovery. Future work considers using the approach in
conjunction with dimensionality reduction techniques to develop parameterized
low-order models for parameterized PDEs. An interesting challenge that we
foresee following that path is that, not only the dynamics may depend on the
parameter, but also the low dimensional embedding.

###### Acknowledgements.

We gratefully acknowledge S. Brunton, U. Fasel, M. Bahamondes, D. Delgado, and
N. Torres for helpful comments and insightful discussions. This work was
funded by U. of Chile internal grant U-Inicia-003/21 and ANID Fondecyt
11220465.

## References

  * Brunton  _et al._ [2016a] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113, 3932 (2016a).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Reviews of modern physics 65, 851 (1993).
  * Benner  _et al._ [2015] P. Benner, S. Gugercin, and K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Review 57, 483 (2015).
  * Rowley and Dawson [2017] C. W. Rowley and S. T. Dawson, Model reduction for flow analysis and control, Annual Review of Fluid Mechanics 49, 387 (2017).
  * Brunton and Noack [2015] S. L. Brunton and B. R. Noack, Closed-loop turbulence control: Progress and challenges, Applied Mechanics Reviews 67, 050801 (2015).
  * Herrmann  _et al._ [2020a] B. Herrmann, M. Behzad, J. M. Cardemil, W. R. Calderón-Muñoz, and R. M. Fernández, Conjugate heat transfer model for feedback control and state estimation in a volumetric solar receiver, Solar Energy 198, 343 (2020a).
  * Herrmann  _et al._ [2022] B. Herrmann, S. L. Brunton, J. E. Pohl, and R. Semaan, Gust mitigation through closed-loop control. ii. feedforward and feedback control, Physical Review Fluids 7, 024706 (2022).
  * Gunzburger [1999] M. D. Gunzburger, Sensitivities, adjoints and flow optimization, International Journal for Numerical Methods in Fluids 31, 53 (1999).
  * Giles and Pierce [2000] M. B. Giles and N. A. Pierce, An introduction to the adjoint approach to design, Flow, Turbulence and Combustion 65, 393 (2000).
  * Herrmann-Priesnitz  _et al._ [2016] B. Herrmann-Priesnitz, W. R. Calderón-Muñoz, A. Valencia, and R. Soto, Thermal design exploration of a swirl flow microchannel heat sink for high heat flux applications based on numerical simulations, Applied Thermal Engineering 109, 22 (2016).
  * Herrmann-Priesnitz  _et al._ [2018] B. Herrmann-Priesnitz, W. R. Calderón-Muñoz, G. Diaz, and R. Soto, Heat transfer enhancement strategies in a swirl flow minichannel heat sink based on hydrodynamic receptivity, International Journal of Heat and Mass Transfer 127, 245 (2018).
  * Grieves and Vickers [2017] M. Grieves and J. Vickers, Digital twin: Mitigating unpredictable, undesirable emergent behavior in complex systems, in  _Transdisciplinary perspectives on complex systems: New findings and approaches_ (Springer, , 2017) pp. 85–113.
  * Hartmann  _et al._ [2018] D. Hartmann, M. Herz, and U. Wever, Model order reduction a key technology for digital twins, in  _Reduced-order modeling (ROM) for simulation and optimization_ (Springer, , 2018) pp. 167–179.
  * Niederer  _et al._ [2021] S. A. Niederer, M. S. Sacks, M. Girolami, and K. Willcox, Scaling digital twins from the artisanal to the industrial, Nature Computational Science 1, 313 (2021).
  * Kapteyn  _et al._ [2021] M. G. Kapteyn, J. V. Pretorius, and K. E. Willcox, A probabilistic graphical model foundation for enabling predictive digital twins at scale, Nature Computational Science 1, 337 (2021).
  * McClellan  _et al._ [2022] A. McClellan, J. Lorenzetti, M. Pavone, and C. Farhat, A physics-based digital twin for model predictive control of autonomous unmanned aerial vehicle landing, Philosophical Transactions of the Royal Society A 380, 20210204 (2022).
  * Ali  _et al._ [2023] W. A. Ali, M. P. Fanti, M. Roccotelli, and L. Ranieri, A review of digital twin technology for electric and autonomous vehicles, Applied Sciences 13, 5871 (2023).
  * Arzani and Dawson [2021] A. Arzani and S. T. Dawson, Data-driven cardiovascular flow modelling: examples and opportunities, Journal of the Royal Society Interface 18, 20200802 (2021).
  * Brunton and Kutz [2019] S. L. Brunton and J. N. Kutz,  _Data-driven science and engineering: Machine learning, dynamical systems, and control_ (Cambridge University Press, , 2019).
  * Gottwald and Reich [2021a] G. A. Gottwald and S. Reich, Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation, Physica D: Nonlinear Phenomena 423, 132911 (2021a).
  * Gottwald and Reich [2021b] G. A. Gottwald and S. Reich, Combining machine learning and data assimilation to forecast dynamical systems from noisy partial observations, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021b).
  * Pathak  _et al._ [2018] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach, Physical review letters 120, 024102 (2018).
  * Vlachas  _et al._ [2018] P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, and P. Koumoutsakos, Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks, Proceedings of the Royal Society A 474, 20170844 (2018).
  * Chen  _et al._ [2018] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, Neural ordinary differential equations, Advances in Neural Information Processing Systems 31 (2018).
  * Raissi  _et al._ [2019] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics 378, 686 (2019).
  * Li  _et al._ [2020] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv preprint arXiv:2010.08895 (2020).
  * Rackauckas  _et al._ [2020] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman, Universal differential equations for scientific machine learning, arXiv preprint arXiv:2001.04385 (2020).
  * Lu  _et al._ [2021] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, Learning nonlinear operators via deeponet based on the universal approximation theorem of operators, Nature machine intelligence 3, 218 (2021).
  * Udrescu and Tegmark [2020] S.-M. Udrescu and M. Tegmark, Ai feynman: A physics-inspired method for symbolic regression, Science Advances 6, eaay2631 (2020).
  * Bongard and Lipson [2007] J. Bongard and H. Lipson, Automated reverse engineering of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 104, 9943 (2007).
  * Schmidt and Lipson [2009] M. Schmidt and H. Lipson, Distilling free-form natural laws from experimental data, science 324, 81 (2009).
  * Cranmer  _et al._ [2020] M. Cranmer, A. Sanchez Gonzalez, P. Battaglia, R. Xu, K. Cranmer, D. Spergel, and S. Ho, in  _Discovering symbolic models from deep learning with inductive biases_ , Vol. 33 (Curran Associates, Inc., 2020) pp. 17429–17442.
  * Raissi  _et al._ [2017] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Machine learning of linear differential equations using gaussian processes, Journal of Computational Physics 348, 683 (2017).
  * Billings [2013] S. A. Billings,  _Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains_ (John Wiley & Sons, , 2013).
  * Schmid [2010] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656, 5 (2010).
  * Rowley  _et al._ [2009] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson, Spectral analysis of nonlinear flows, Journal of Fluid Mechanics 641, 115 (2009).
  * Kutz  _et al._ [2016] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor,  _Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems_ (SIAM, , 2016).
  * Schmid [2022] P. J. Schmid, Dynamic mode decomposition and its variants, Annual Review of Fluid Mechanics 54, 225 (2022).
  * Herrmann  _et al._ [2021] B. Herrmann, P. J. Baddoo, R. Semaan, S. L. Brunton, and B. J. McKeon, Data-driven resolvent analysis, Journal of Fluid Mechanics 918, A10 (2021).
  * Baddoo  _et al._ [2022] P. J. Baddoo, B. Herrmann, B. J. McKeon, and S. L. Brunton, Kernel learning for robust dynamic mode decomposition: linear and nonlinear disambiguation optimization, Proceedings of the Royal Society A 478, 20210830 (2022).
  * Baddoo  _et al._ [2023] P. J. Baddoo, B. Herrmann, B. J. McKeon, J. Nathan Kutz, and S. L. Brunton, Physics-informed dynamic mode decomposition, Proceedings of the Royal Society A 479, 20220576 (2023).
  * Mezić [2005] I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics 41, 309 (2005).
  * Mezić [2013] I. Mezić, Analysis of fluid flows via spectral properties of the koopman operator, Annual Review of Fluid Mechanics 45, 357 (2013).
  * Brunton  _et al._ [2022] S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz, Modern koopman theory for dynamical systems, SIAM Review 64, 229 (2022).
  * Peherstorfer and Willcox [2016] B. Peherstorfer and K. Willcox, Data-driven operator inference for nonintrusive projection-based model reduction, Computer Methods in Applied Mechanics and Engineering 306, 196 (2016).
  * Qian  _et al._ [2020] E. Qian, B. Kramer, B. Peherstorfer, and K. Willcox, Lift & learn: Physics-informed machine learning for large-scale nonlinear dynamical systems, Physica D: Nonlinear Phenomena 406, 132401 (2020).
  * Cenedese  _et al._ [2022a] M. Cenedese, J. Axås, B. Bäuerlein, K. Avila, and G. Haller, Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds, Nature communications 13, 872 (2022a).
  * Cenedese  _et al._ [2022b] M. Cenedese, J. Axås, H. Yang, M. Eriten, and G. Haller, Data-driven nonlinear model reduction to spectral submanifolds in mechanical systems, Philosophical Transactions of the Royal Society A 380, 20210194 (2022b).
  * Axås  _et al._ [2023] J. Axås, M. Cenedese, and G. Haller, Fast data-driven model reduction for nonlinear dynamical systems, Nonlinear Dynamics 111, 7941 (2023).
  * Rudy  _et al._ [2017] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Data-driven discovery of partial differential equations, Science Advances 3, e1602614 (2017).
  * Tibshirani [1996] R. Tibshirani, Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society Series B: Statistical Methodology 58, 267 (1996).
  * Joshi  _et al._ [2022] C. Joshi, S. Ray, L. M. Lemma, M. Varghese, G. Sharp, Z. Dogic, A. Baskaran, and M. F. Hagan, Data-driven discovery of active nematic hydrodynamics, Physical review letters 129, 258001 (2022).
  * Golden  _et al._ [2023] M. Golden, R. O. Grigoriev, J. Nambisan, and A. Fernandez-Nieves, Physically informed data-driven modeling of active nematics, Science Advances 9, eabq6120 (2023).
  * Rubio-Herrero  _et al._ [2022] J. Rubio-Herrero, C. O. Marrero, and W.-T. L. Fan, Modeling atmospheric data and identifying dynamics temporal data-driven modeling of air pollutants, Journal of Cleaner Production 333, 129863 (2022).
  * Pasquato  _et al._ [2022] M. Pasquato, M. Abbas, A. A. Trani, M. Nori, J. A. Kwiecinski, P. Trevisan, V. F. Braga, G. Bono, and A. V. Macciò, Sparse identification of variable star dynamics, The Astrophysical Journal 930, 161 (2022).
  * Mangan  _et al._ [2016] N. M. Mangan, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Inferring biological networks by sparse identification of nonlinear dynamics, IEEE Transactions on Molecular, Biological and Multi-Scale Communications 2, 52 (2016).
  * Lagergren  _et al._ [2020] J. H. Lagergren, J. T. Nardini, G. Michael Lavigne, E. M. Rutter, and K. B. Flores, Learning partial differential equations for biological transport models from noisy spatio-temporal data, Proceedings of the Royal Society A 476, 20190800 (2020).
  * Boninsegna  _et al._ [2018] L. Boninsegna, F. Nüske, and C. Clementi, Sparse learning of stochastic dynamical equations, The Journal of chemical physics 148 (2018).
  * Hoffmann  _et al._ [2019] M. Hoffmann, C. Fröhner, and F. Noé, Reactive sindy: Discovering governing reactions from concentration data, The Journal of Chemical Physics 150 (2019).
  * Bhadriraju  _et al._ [2020] B. Bhadriraju, M. S. F. Bangi, A. Narasingam, and J. S.-I. Kwon, Operable adaptive sparse identification of systems: Application to chemical processes, AIChE Journal 66, e16980 (2020).
  * Scheffold  _et al._ [2021] L. Scheffold, T. Finkler, and U. Piechottka, Gray-box system modeling using symbolic regression and nonlinear model predictive control of a semibatch polymerization, Computers & Chemical Engineering 146, 107204 (2021).
  * Dale and Bhat [2018] R. Dale and H. S. Bhat, Equations of mind: Data science for inferring nonlinear dynamics of socio-cognitive systems, Cognitive Systems Research 52, 275 (2018).
  * Kaheman  _et al._ [2019] K. Kaheman, E. Kaiser, B. Strom, J. N. Kutz, and S. L. Brunton, Learning discrepancy models from experimental data, arXiv preprint arXiv:1909.08574 (2019).
  * de Silva  _et al._ [2020a] B. M. de Silva, D. M. Higdon, S. L. Brunton, and J. N. Kutz, Discovery of physics from data: Universal laws and discrepancies, Frontiers in artificial intelligence 3, 25 (2020a).
  * Ebers  _et al._ [2022] M. R. Ebers, K. M. Steele, and J. N. Kutz, Discrepancy modeling framework: Learning missing physics, modeling systematic residuals, and disambiguating between deterministic and random effects, arXiv preprint arXiv:2203.05164 (2022).
  * Jiang  _et al._ [2021] Y.-X. Jiang, X. Xiong, S. Zhang, J.-X. Wang, J.-C. Li, and L. Du, Modeling and prediction of the transmission dynamics of covid-19 based on the sindy-lm method, Nonlinear Dynamics 105, 2775 (2021).
  * Loiseau and Brunton [2018] J.-C. Loiseau and S. L. Brunton, Constrained sparse galerkin regression, Journal of Fluid Mechanics 838, 42 (2018).
  * Loiseau  _et al._ [2018] J.-C. Loiseau, B. R. Noack, and S. L. Brunton, Sparse reduced-order modelling: sensor-based dynamics to full-state estimation, Journal of Fluid Mechanics 844, 459 (2018).
  * El Sayed M  _et al._ [2018] Y. El Sayed M, R. Semaan, and R. Radespiel, Sparse modeling of the lift gains of a high-lift configuration with periodic coanda blowing, in  _2018 AIAA Aerospace Sciences Meeting_ (2018) p. 1054.
  * Chang and Zhang [2019] H. Chang and D. Zhang, Machine learning subsurface flow equations from data, Computational Geosciences 23, 895 (2019).
  * Deng  _et al._ [2020] N. Deng, B. R. Noack, M. Morzyński, and L. R. Pastur, Low-order model for successive bifurcations of the fluidic pinball, Journal of fluid mechanics 884, A37 (2020).
  * Loiseau [2020] J.-C. Loiseau, Data-driven modeling of the chaotic thermal convection in an annular thermosyphon, Theoretical and Computational Fluid Dynamics 34, 339 (2020).
  * Fukami  _et al._ [2021] K. Fukami, T. Murata, K. Zhang, and K. Fukagata, Sparse identification of nonlinear dynamics with low-dimensionalized flow representations, Journal of Fluid Mechanics 926, A10 (2021).
  * Guan  _et al._ [2021] Y. Guan, S. L. Brunton, and I. Novosselov, Sparse nonlinear models of chaotic electroconvection, Royal Society Open Science 8, 202367 (2021).
  * Deng  _et al._ [2021] N. Deng, B. R. Noack, M. Morzyński, and L. R. Pastur, Galerkin force model for transient and post-transient dynamics of the fluidic pinball, Journal of Fluid Mechanics 918, A4 (2021).
  * Callaham  _et al._ [2022a] J. L. Callaham, G. Rigas, J.-C. Loiseau, and S. L. Brunton, An empirical mean-field model of symmetry-breaking in a turbulent wake, Science Advances 8, eabm4786 (2022a).
  * Callaham  _et al._ [2022b] J. L. Callaham, S. L. Brunton, and J.-C. Loiseau, On the role of nonlinear correlations in reduced-order modelling, Journal of Fluid Mechanics 938, A1 (2022b).
  * Khoo  _et al._ [2022] Z. C. Khoo, C. H. Chan, and Y. Hwang, A sparse optimal closure for a reduced-order model of wall-bounded turbulence, Journal of Fluid Mechanics 939, A11 (2022).
  * Sansica  _et al._ [2022] A. Sansica, J.-C. Loiseau, M. Kanamori, A. Hashimoto, and J.-C. Robinet, System identification of two-dimensional transonic buffet, AIAA Journal 60, 3090 (2022).
  * Xiao  _et al._ [2023] Q. Xiao, J. Wang, X. Yang, and B. Jiang, Construction of a reduced-order model of an electro-osmotic micromixer and discovery of attractors for petal structure, Physics of Fluids 35 (2023).
  * Subramanian  _et al._ [2021] R. Subramanian, R. R. Moar, and S. Singh, White-box machine learning approaches to identify governing equations for overall dynamics of manufacturing systems: A case study on distillation column, Machine Learning with Applications 3, 100014 (2021).
  * Ito  _et al._ [2023] Y. Ito, T. Hato, and A. Kano, Physics-informed machine learning for surrogate modeling of heat transfer phenomena, Journal of Computational and Nonlinear Dynamics 18, 111001 (2023).
  * Thaler  _et al._ [2019] S. Thaler, L. Paehler, and N. A. Adams, Sparse identification of truncation errors, Journal of Computational Physics 397, 108851 (2019).
  * Sorokina  _et al._ [2016] M. Sorokina, S. Sygletos, and S. Turitsyn, Sparse identification for nonlinear optical communication systems: Sino method, Optics Express 24, 30433 (2016).
  * Dam  _et al._ [2017] M. Dam, M. Brøns, J. Juul Rasmussen, V. Naulin, and J. S. Hesthaven, Sparse identification of a predator-prey system from simulation data of a convection model, Physics of Plasmas 24 (2017).
  * Kaptanoglu  _et al._ [2021a] A. A. Kaptanoglu, K. D. Morgan, C. J. Hansen, and S. L. Brunton, Physics-constrained, low-dimensional models for magnetohydrodynamics: First-principles and data-driven approaches, Physical Review E 104, 015206 (2021a).
  * Alves and Fiuza [2022] E. P. Alves and F. Fiuza, Data-driven discovery of reduced plasma physics models from fully kinetic simulations, Physical Review Research 4, 033192 (2022).
  * Stanković  _et al._ [2020] A. M. Stanković, A. A. Sarić, A. T. Sarić, and M. K. Transtrum, Data-driven symbolic regression for identification of nonlinear dynamics in power systems, in  _2020 IEEE Power & Energy Society General Meeting_ (IEEE, 2020) pp. 1–5.
  * Cai  _et al._ [2022] Y. Cai, X. Wang, G. Joós, and I. Kamwa, An online data-driven method to locate forced oscillation sources from power plants based on sparse identification of nonlinear dynamics (sindy), IEEE Transactions on Power Systems 38, 2085 (2022).
  * Lai and Nagarajaiah [2019] Z. Lai and S. Nagarajaiah, Sparse structural system identification method for nonlinear dynamic systems with hysteresis/inelastic behavior, Mechanical Systems and Signal Processing 117, 813 (2019).
  * Chatterjee  _et al._ [2023] T. Chatterjee, A. D. Shaw, M. I. Friswell, and H. H. Khodaparast, Sparse bayesian machine learning for the interpretable identification of nonlinear structural dynamics: Towards the experimental data-driven discovery of a quasi zero stiffness device, Mechanical Systems and Signal Processing 205, 110858 (2023).
  * Beetham and Capecelatro [2020] S. Beetham and J. Capecelatro, Formulating turbulence closures using sparse regression with embedded form invariance, Physical Review Fluids 5, 084611 (2020).
  * Schmelzer  _et al._ [2020] M. Schmelzer, R. P. Dwight, and P. Cinnella, Discovery of algebraic reynolds-stress models using sparse symbolic regression, Flow, Turbulence and Combustion 104, 579 (2020).
  * Zanna and Bolton [2020] L. Zanna and T. Bolton, Data-driven equation discovery of ocean mesoscale closures, Geophysical Research Letters 47, e2020GL088376 (2020).
  * Beetham  _et al._ [2021] S. Beetham, R. O. Fox, and J. Capecelatro, Sparse identification of multiphase turbulence closures for coupled fluid–particle flows, Journal of Fluid Mechanics 914, A11 (2021).
  * Schaeffer [2017] H. Schaeffer, Learning partial differential equations via data discovery and sparse optimization, Proceedings of the Royal Society A 473, 20160446 (2017).
  * Mangan  _et al._ [2019] N. M. Mangan, T. Askham, S. L. Brunton, J. N. Kutz, and J. L. Proctor, Model selection for hybrid dynamical systems via sparse regression, Proceedings of the Royal Society A 475, 20180534 (2019).
  * Klimovskaia  _et al._ [2016] A. Klimovskaia, S. Ganscha, and M. Claassen, Sparse regression based structure learning of stochastic reaction networks from single cell snapshot time series, PLoS computational biology 12, e1005234 (2016).
  * Brückner  _et al._ [2020] D. B. Brückner, P. Ronceray, and C. P. Broedersz, Inferring the dynamics of underdamped stochastic systems, Physical review letters 125, 058103 (2020).
  * Dai  _et al._ [2020] M. Dai, T. Gao, Y. Lu, Y. Zheng, and J. Duan, Detecting the maximum likelihood transition path from data of stochastic dynamical systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (2020).
  * Callaham  _et al._ [2021] J. L. Callaham, J.-C. Loiseau, G. Rigas, and S. L. Brunton, Nonlinear stochastic modelling with langevin regression, Proceedings of the Royal Society A 477, 20210092 (2021).
  * Huang  _et al._ [2022] Y. Huang, Y. Mabrouk, G. Gompper, and B. Sabass, Sparse inference and active learning of stochastic differential equations from data, Scientific Reports 12, 21691 (2022).
  * Brunton  _et al._ [2016b] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Sparse identification of nonlinear dynamics with control (sindyc), IFAC-PapersOnLine 49, 710 (2016b).
  * Kaiser  _et al._ [2018] E. Kaiser, J. N. Kutz, and S. L. Brunton, Sparse identification of nonlinear dynamics for model predictive control in the low-data limit, Proceedings of the Royal Society A 474, 20180335 (2018).
  * Fasel  _et al._ [2021] U. Fasel, E. Kaiser, J. N. Kutz, B. W. Brunton, and S. L. Brunton, Sindy with control: A tutorial, in  _2021 60th IEEE Conference on Decision and Control (CDC)_ (IEEE, 2021) pp. 16–21.
  * Kaheman  _et al._ [2020] K. Kaheman, J. N. Kutz, and S. L. Brunton, Sindy-pi: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics, Proceedings of the Royal Society A 476, 20200279 (2020).
  * Kaptanoglu  _et al._ [2021b] A. A. Kaptanoglu, J. L. Callaham, A. Aravkin, C. J. Hansen, and S. L. Brunton, Promoting global stability in data-driven models of quadratic nonlinear dynamics, Physical Review Fluids 6, 094401 (2021b).
  * Bakarji  _et al._ [2022] J. Bakarji, J. Callaham, S. L. Brunton, and J. N. Kutz, Dimensionally consistent learning with buckingham pi, Nature Computational Science 2, 834 (2022).
  * Zhang and Lin [2018] S. Zhang and G. Lin, Robust data-driven discovery of governing physical laws with error bars, Proceedings of the Royal Society A 474, 20180305 (2018).
  * Ram  _et al._ [2021] P. R. M. Ram, U. Römer, and R. Semaan, Bayesian dynamical system identification with unified sparsity priors and model uncertainty, arXiv preprint arXiv:2103.05090 (2021).
  * Hirsh  _et al._ [2022] S. M. Hirsh, D. A. Barajas-Solano, and J. N. Kutz, Sparsifying priors for bayesian uncertainty quantification in model discovery, Royal Society Open Science 9, 211823 (2022).
  * Fasel  _et al._ [2022] U. Fasel, J. N. Kutz, B. W. Brunton, and S. L. Brunton, Ensemble-sindy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control, Proceedings of the Royal Society A 478, 20210904 (2022).
  * Gao  _et al._ [2023] L. Gao, U. Fasel, S. L. Brunton, and J. N. Kutz, Convergence of uncertainty estimates in ensemble and bayesian sparse model discovery, arXiv preprint arXiv:2301.12649 (2023).
  * Quade  _et al._ [2018] M. Quade, M. Abel, J. Nathan Kutz, and S. L. Brunton, Sparse identification of nonlinear dynamics for rapid model recovery, Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (2018).
  * Schaeffer  _et al._ [2018] H. Schaeffer, G. Tran, and R. Ward, Extracting sparse high-dimensional dynamics from limited data, SIAM Journal on Applied Mathematics 78, 3279 (2018).
  * Champion  _et al._ [2019a] K. P. Champion, S. L. Brunton, and J. N. Kutz, Discovery of nonlinear multiscale systems: Sampling strategies and embeddings, SIAM Journal on Applied Dynamical Systems 18, 312 (2019a).
  * Mangan  _et al._ [2017] N. M. Mangan, J. N. Kutz, S. L. Brunton, and J. L. Proctor, Model selection for dynamical systems via sparse regression and information criteria, Proceedings of the Royal Society A 473, 20170009 (2017).
  * Dong  _et al._ [2023] X. Dong, Y.-L. Bai, Y. Lu, and M. Fan, An improved sparse identification of nonlinear dynamics with akaike information criterion and group sparsity, Nonlinear Dynamics 111, 1485 (2023).
  * Schaeffer and McCalla [2017] H. Schaeffer and S. G. McCalla, Sparse model selection via integral terms, Physical Review E 96, 023302 (2017).
  * Gurevich  _et al._ [2019] D. R. Gurevich, P. A. Reinbold, and R. O. Grigoriev, Robust and optimal sparse regression for nonlinear pde models, Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (2019).
  * Messenger and Bortz [2021a] D. A. Messenger and D. M. Bortz, Weak sindy: Galerkin-based data-driven model selection, Multiscale Modeling & Simulation 19, 1474 (2021a).
  * Messenger and Bortz [2021b] D. A. Messenger and D. M. Bortz, Weak sindy for partial differential equations, Journal of Computational Physics 443, 110525 (2021b).
  * Delahunt and Kutz [2022] C. B. Delahunt and J. N. Kutz, A toolkit for data-driven discovery of governing equations in high-noise regimes, IEEE Access 10, 31210 (2022).
  * Kaheman  _et al._ [2022] K. Kaheman, S. L. Brunton, and J. N. Kutz, Automatic differentiation to simultaneously identify nonlinear dynamics and extract noise probability distributions from data, Machine Learning: Science and Technology 3, 015031 (2022).
  * Tran and Ward [2017] G. Tran and R. Ward, Exact recovery of chaotic systems from highly corrupted data, Multiscale Modeling & Simulation 15, 1108 (2017).
  * Champion  _et al._ [2020] K. Champion, P. Zheng, A. Y. Aravkin, S. L. Brunton, and J. N. Kutz, A unified sparse optimization framework to learn parsimonious physics-informed models from data, IEEE Access 8, 169259 (2020).
  * Reinbold  _et al._ [2020] P. A. Reinbold, D. R. Gurevich, and R. O. Grigoriev, Using noisy or incomplete data to discover models of spatiotemporal dynamics, Physical Review E 101, 010203 (2020).
  * Reinbold  _et al._ [2021] P. A. Reinbold, L. M. Kageorge, M. F. Schatz, and R. O. Grigoriev, Robust learning from noisy, incomplete, high-dimensional experimental data via physically constrained symbolic regression, Nature communications 12, 3219 (2021).
  * Somacal  _et al._ [2022] A. Somacal, Y. Barrera, L. Boechi, M. Jonckheere, V. Lefieux, D. Picard, and E. Smucler, Uncovering differential equations from data with hidden variables, Physical Review E 105, 054209 (2022).
  * Bakarji  _et al._ [2023] J. Bakarji, K. Champion, J. Nathan Kutz, and S. L. Brunton, Discovering governing equations from partial measurements with deep delay autoencoders, Proceedings of the Royal Society A 479, 20230422 (2023).
  * de Silva  _et al._ [2020b] B. M. de Silva, K. Champion, M. Quade, J.-C. Loiseau, J. N. Kutz, and S. L. Brunton, Pysindy: A python package for the sparse identification of nonlinear dynamical systems from data, Journal of Open Source Software 5, 2104 (2020b).
  * Kaptanoglu  _et al._ [2022] A. A. Kaptanoglu, B. M. de Silva, U. Fasel, K. Kaheman, A. J. Goldschmidt, J. L. Callaham, C. B. Delahunt, Z. G. Nicolaou, K. Champion, J.-C. Loiseau, _et al._ , Pysindy: A comprehensive python package for robust sparse system identification, Journal of Open Source Software 7, 3994 (2022).
  * Kaptanoglu  _et al._ [2023] A. A. Kaptanoglu, L. Zhang, Z. G. Nicolaou, U. Fasel, and S. L. Brunton, Benchmarking sparse system identification with low-dimensional chaos, Nonlinear Dynamics , 1 (2023).
  * Gilpin [2021] W. Gilpin, Chaos as an interpretable benchmark for forecasting and data-driven modelling, Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks 1 (2021).
  * Champion  _et al._ [2019b] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton, Data-driven discovery of coordinates and governing equations, Proceedings of the National Academy of Sciences 116, 22445 (2019b).
  * Fries  _et al._ [2022] W. D. Fries, X. He, and Y. Choi, Lasdi: Parametric latent space dynamics identification, Computer Methods in Applied Mechanics and Engineering 399, 115436 (2022).
  * Conti  _et al._ [2023] P. Conti, G. Gobat, S. Fresca, A. Manzoni, and A. Frangi, Reduced order modeling of parametrized systems through autoencoders and sindy approach: continuation of periodic solutions, Computer Methods in Applied Mechanics and Engineering 411, 116072 (2023).
  * Meidani and Farimani [2023] K. Meidani and A. B. Farimani, Identification of parametric dynamical systems using integer programming, Expert Systems with Applications 219, 119622 (2023).
  * Schaeffer  _et al._ [2017] H. Schaeffer, G. Tran, and R. Ward, Learning dynamical systems and bifurcation via group sparsity, arXiv preprint arXiv:1709.01558 (2017).
  * Rudy  _et al._ [2019] S. Rudy, A. Alla, S. L. Brunton, and J. N. Kutz, Data-driven identification of parametric partial differential equations, SIAM Journal on Applied Dynamical Systems 18, 643 (2019).
  * Nicolaou  _et al._ [2023] Z. G. Nicolaou, G. Huo, Y. Chen, S. L. Brunton, and J. N. Kutz, Data-driven discovery and extrapolation of parameterized pattern-forming dynamics, Physical Review Research 5, L042017 (2023).
  * Dankowicz and Schilder [2013] H. Dankowicz and F. Schilder,  _Recipes for continuation_ (SIAM, , 2013).
  * Ludwig  _et al._ [1978] D. Ludwig, D. D. Jones, C. S. Holling, _et al._ , Qualitative analysis of insect outbreak systems: the spruce budworm and forest, Journal of animal ecology 47, 315 (1978).
  * Herrmann  _et al._ [2020b] B. Herrmann, P. Oswald, R. Semaan, and S. L. Brunton, Modeling synchronization in forced turbulent oscillator flows, Communications Physics 3, 1 (2020b).
  * Lorenz [1963] E. N. Lorenz, Deterministic nonperiodic flow, Journal of atmospheric sciences 20, 130 (1963).