SINDy-PI: A Robust Algorithm for Parallel Implicit Sparse Identification of Nonlinear Dynamics

04/05/2020
by   Kadierdan Kaheman, et al.
University of Washington
5

Accurately modeling the nonlinear dynamics of a system from measurement data is a challenging yet vital topic. The sparse identification of nonlinear dynamics (SINDy) algorithm is one approach to discover dynamical systems models from data. Although extensions have been developed to identify implicit dynamics, or dynamics described by rational functions, these extensions are extremely sensitive to noise. In this work, we develop SINDy-PI (parallel, implicit), a robust variant of the SINDy algorithm to identify implicit dynamics and rational nonlinearities. The SINDy-PI framework includes multiple optimization algorithms and a principled approach to model selection. We demonstrate the ability of this algorithm to learn implicit ordinary and partial differential equations and conservation laws from limited and noisy data. In particular, we show that the proposed approach is several orders of magnitude more noise robust than previous approaches, and may be used to identify a class of complex ODE and PDE dynamics that were previously unattainable with SINDy, including for the double pendulum dynamics and the Belousov Zhabotinsky (BZ) reaction.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 7

page 12

11/22/2021

Ensemble-SINDy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control

Sparse model identification enables the discovery of nonlinear dynamical...
11/12/2021

PySINDy: A comprehensive Python package for robust sparse system identification

Automated data-driven modeling, the process of directly discovering the ...
09/11/2021

Structure-preserving Sparse Identification of Nonlinear Dynamics for Data-driven Modeling

Discovery of dynamical systems from data forms the foundation for data-d...
01/07/2021

CINDy: Conditional gradient-based Identification of Non-linear Dynamics – Noise-robust recovery

Governing equations are essential to the study of nonlinear dynamics, of...
09/18/2019

Learning Discrepancy Models From Experimental Data

First principles modeling of physical systems has led to significant tec...
07/17/2019

Stability selection enables robust learning of partial differential equations from limited noisy data

We present a statistical learning framework for robust identification of...
05/09/2020

Weak SINDy: A Data-Driven Galerkin Method for System Identification

We present a weak formulation and discretization of the system discovery...

Code Repositories

SINDy-PI

SINDy-PI: A Robust Algorithm for Parallel Implicit Sparse Identification of Nonlinear Dynamics


view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Discovering dynamical system models from data is critically important across science and engineering. Traditionally, models are derived from first principles, although this approach may be prohibitively challenging in many fields, such as climate science, finance, and biology. Fortunately, data-driven model discovery (i.e., system identification) is a rapidly developing field [Brunton2019book], with a range of techniques including classic linear approaches [Nelles2013book, ljung2010arc], dynamic mode decomposition (DMD) [Schmid2010jfm, Kutz2016book] and Koopman theory [Budivsic2012chaos, Mezic2013arfm, Williams2015jnls, klus2017data]

, nonlinear autoregressive models 

[Akaike1969annals, Billings2013book]

, neural networks 

[yang2018physics, Wehmeyer2018jcp, Mardt2018natcomm, vlachas2018data, pathak2018model, lu2019deepxde, Raissi2019jcp, Champion2019pnas, raissi2020science], Gaussian process regression [Raissi2017arxiva, Raissi2017arxiv], nonlinear Laplacian spectral analysis [Giannakis2012pnas], diffusion maps [Yair2017pnas]

, genetic programming 

[Bongard2007pnas, Schmidt2009science, Daniels2015naturecomm], and sparse regression [Brunton2016SINDy, Rudy2017PDE-find], to highlight some of the recent developments. Of particular note is a recent push towards parsimonious modeling [Bongard2007pnas, Schmidt2009science, Brunton2016SINDy], which favors Pareto-optimal models with the lowest complexity required to describe the observed data. These models benefit from being interpretable, and they tend to generalize and prevent overfitting. The sparse identification of nonlinear dynamics (SINDy) algorithm [Brunton2016SINDy] discovers parsimonious models through a sparsity-promoting optimization to select only a few model terms from a library of candidate functions. SINDy has been widely adopted in the community [Sorokina2016oe, Dam2017pf, Schaeffer2017prsa, narasingam2018data, JC2018ConGalerkin, Quade2018chaos, boninsegna2018sparse, Loiseau2018jfm, Zhang2018arxiv, Hoffmann2018arxiv, mangan2019model, Thaler2019jcp, lai2019sparse, de2019discovery, Gelss2019mindy, Guang2018], but it relies on the dynamics having a sparse representation in a pre-defined library, making it difficult to discover implicit dynamics and rational functions. The implicit-SINDy extension [Mangan2016ImSINDY] makes it possible to identify these implicit functions, although this algorithm is extremely sensitive to noise. In this work, we develop a robust, parallel algorithm for the sparse identification of implicit dynamics, making it possible to explore entirely new classes of systems that were previously inaccessible.

Parsimonious modeling has a rich history, with many scientific advances being argued on the basis of Occam’s razor, that the simplest model is likely the correct one. SINDy exemplifies this principle, identifying a potentially nonlinear model with the fewest terms required to describe how the measurement data changes in time. The basic idea behind SINDy may be illustrated on a one-dimensional system ; the general formulation for multidimensional dynamics will be described in the following sections. An interpretable form of the nonlinear dynamics may be learned by writing the rate of change of the state of the system as a sparse linear combination of a few terms in a library of candidate functions, :

(1)

where each is prescribed candidate term (e.g.

). The derivative of the state and the library of candidate functions may both be computed from measured trajectory data. It then remains to solve for a sparse vector

with nonzero entries indicating which functions are active in characterizing the dynamics. The resulting models strike a balance between accuracy and efficiency, and they are highly interpretable by construction. In a short time, the SINDy algorithm has been extended to include inputs and control [Kaiser2018prsa], to identify partial differential equations [Rudy2017PDE-find, Schaeffer2017prsa], to incorporate physically relevant constraints [JC2018ConGalerkin]

, to include tensor bases 

[Gelss2019mindy], and to incorporate integral terms for denoising [Schaeffer2017IntegralSINDy, Reinbold2020pre]. These extensions and its simple formulation in terms of a generalized linear model in (1) have resulted in SINDy being adopted in the fields of fluid mechanics [JC2018ConGalerkin, Loiseau2018jfm], nonlinear optics [Sorokina2016oe], plasma physics [Dam2017pf], chemical reactions [Hoffmann2018arxiv, Boninsegna2018SSR, narasingam2018data], numerical methods [Thaler2019jcp], and structural modeling [lai2019sparse].

The generalized linear model in (1) does not readily lend itself to representing implicit dynamics and rational functions, which are not naturally expressible as sum of a few basis functions. Instead, the implicit-SINDy algorithm [Mangan2016ImSINDY] reformulates the SINDy problem in an implicit form:

(2)

This formulation is flexible enough to handle a much broader class of dynamics with rational function nonlinearities, such as which may be rewritten as . However, the sparsest vector that satisfies (2) is the trivial solution . Thus, the implicit-SINDy algorithm leverages a recent non-convex optimization procedure [Wright2009ieeetpami, Qu2014ADM] to find the sparsest vector in the null space of . For even small amounts of noise, the dimension of the null space will become prohibitively large, making this approach extremely sensitive to noise and compromising the model discovery process.

This work develops an optimization and model selection framework that recasts implicit-SINDy as a convex problem, making it as noise robust as the original non-implicit SINDy algorithm and enabling the identification of implicit ODEs and PDEs that were previously inaccessible. The key to making the implicit-SINDy algorithm robust is the realization that if we know even a single term in the dynamics, corresponding to a non-zero entry , then we can rewrite (2) in a non-implicit form

(3)

where and have the -th element removed. Because none of these terms are known a priori, we sweep through the library, term by term, testing (3) for a sparse model that fits the data. This procedure is highly parallelizable and provides critical information for model selection. Our approach is related to the recent work of Zhang et al. [Guang2018], which also makes the implicit problem more robust by testing candidate functions individually. However, there are a number of key differences in the present approach. Our work explicitly considers rational nonlinearities to discover exceedingly complex implicit PDEs, such as the Belousov-Zhabotinsky (BZ) reaction. Our framework also provides several new greedy algorithms, including parallel and constrained formulations. We further extend this method to include the effect of control inputs, making it applicable to robotic systems [Koryakovskiy2018], and we use this procedure to discover complex Hamiltonians. Finally, our approach provides guidance on model selection, a comprehensive comparison with previous methods, and a careful analysis of noise robustness.

2 Background

We briefly introduce the full multidimensional SINDy and implicit-SINDy algorithms, which will provide a foundation for our robust implicit identification algorithm in Sec. 3.

2.1 Sparse Identification of Nonlinear Dynamics

The goal of SINDy [Brunton2016SINDy] is to discover a dynamical system

(4)

from time-series data of the state . We assume that the dynamics, encoded by the function , admit a sparse representation in a library of candidate functions:

(5)

Thus, each row equation in (4) may be written as

(6)

where is a sparse vector, indicating which terms are active in the dynamics.

We determine the nonzero entries of through sparse regression based on trajectory data. The time-series data is arranged into a matrix , and the associated time derivative matrix is computed using an appropriate numerical differentiation scheme [Chartrand2011TVReg, Brunton2016SINDy, Brunton2019book]. It is then possible to evaluate the library on trajectory data in so that each column of is a function evaluated on the snapshots in .

It is now possible to write the dynamical system in terms of a generalized linear model, evaluated on trajectory data:

(7)

There are several approaches to identify the sparse matrix of coefficients , including sequentially thresholded least squares (STLSQ) [Brunton2016SINDy, zhang2019convergence], LASSO [Tibshirani1994], sparse relaxed regularized regression (SR3) [ZhengSR3, Champion2019SINDySR3], stepwise sparse regression (SSR) [Boninsegna2018SSR], and Bayesian approaches [Guang2018, Pan2016BayesianSINDy]. It is possible to augment the library to include partial derivatives for the identification of partial differential equations (PDEs) [Rudy2017PDE-find, Schaeffer2017prsa]. Similarly, it is possible to include external forcing terms in the library , enabling the identification of forced and actively controlled systems [Kaiser2018prsa]. To alleviate the effect of noise, it is possible to reframe the SINDy problem in terms of an integral formulation [Schaeffer2017IntegralSINDy, Reinbold2020pre].

2.2 Implicit Sparse Identification of Nonlinear Dynamics

The implicit-SINDy algorithm [Mangan2016ImSINDY] extends SINDy to identify implicit differential equations

(8)

and in particular, systems that include rational functions in the dynamics, such as chemical reactions and metabolic networks that have a separation of timescales.

The implicit-SINDy generalizes the library in (7) to include functions of and :

(9)

However, this approach requires solving for a matrix whose columns are sparse vectors in the null space of . This approach is non-convex, relying on the alternating directions method (ADM) [Mangan2016ImSINDY, Qu2014ADM], and null space computations are highly ill-conditioned for noisy data [Gavish2014ieeetit, Mangan2016ImSINDY, Brunton2019book], thus inspiring the current work and mathematical innovations.

Figure 1: The illustration of the SINDy-PI algorithm on Michaelis-Menten dynamics. (a) The Michaelis-Menten system is simulated, and measurement data is provided to SINDy-PI. (b) Multiple possible left-hand side functions are tested at the same time. (c) The candidate model prediction error is calculated, and the best model is selected.

3 SINDy-PI: Robust Parallel Identification of Implicit Dynamics

We have developed the SINDy-PI (parallel, implicit) framework for the robust identification of implicit dynamics, bypassing the null space approach discussed in Sec. 2.2. The idea is that if even a single term in the dynamics (8) is known, it is possible to rewrite (9) as

(10)

where is the library with the column removed. Equation (10) is no longer in implicit form, and the sparse coefficient matrix corresponding to the remaining terms may be solved for using previously developed SINDy techniques [Brunton2016SINDy, ZhengSR3, Champion2019SINDySR3, Boninsegna2018SSR, Guang2018, Pan2016BayesianSINDy, Rudy2017PDE-find, Schaeffer2017prsa, Schaeffer2017IntegralSINDy, Reinbold2020pre]. Because there is no null space calculation, the resulting algorithm is considerably more robust to noise than the implicit-SINDy algorithm [Mangan2016ImSINDY], i.e. we longer have to deal with an ill-conditioned null space problem. However, in general, the entire point of SINDy is that the dynamics are not known ahead of time, and so it is necessary to test each candidate function until one of the models in (10) admits a sparse and accurate solution. When an incorrect candidate term is used, then the algorithm results in a dense (non-sparse) model and an inaccurate model fit, and when a correct term is included, the algorithm identifies a sparse model and an accurate model fit. In this way, it is clear when the algorithm has identified the correct model. Moreover, there is a wealth of redundant information, since each term in the correct model may be used as the candidate function on the left hand side, and the resulting models may be cross-referenced. This approach is highly parallelizable, and each candidate term may be tested simultaneously in parallel. The non-parallel formulation in (10) was recently introduced by Zhang et al. [Guang2018] in the context of Bayesian regression, where they also make the implicit problem more robust by testing candidate functions individually; however, they do not consider dynamics with rational function nonlinearities or control inputs. In this work, we extend the robust implicit formulation to identify several challenging implicit ODE and PDE systems with rational function nonlinearities, which are ubiquitous in engineering and natural systems, and systems with external forcing and control inputs. We also introduce the parallel formulation and model selection frameworks. Further, we will introduce a constrained optimization framework to simultaneously test all candidate functions.

3.1 Model Selection

For each candidate function in (9), we obtain one candidate model. When the candidate function is not in the true dynamics, then the resulting coefficient vector will not be sparse and there will be large prediction error. In contrast, when a correct candidate function is selected, then we obtain a sparse coefficient vector and small prediction error. For an implicit dynamical system, there may be several different implicit equations that must be identified, resulting in several candidate functions that admit sparse models. The sequentially thresholded least squares (STLSQ) algorithm that we use here, and whose convergence properties are considered by Zhang and Schaeffer [zhang2019convergence], iteratively computes a least-squares solution to minimize and then zero out small entries in that are below a set threshold . This threshold

is a hyperparameter that must be tuned to select the model that most accurately balances accuracy and efficiency. Thus, we must employ model selection techniques to identify the implicit models that best supports the data, while remaining as simple as possible.

There are several valid approaches to model selection. To select a parsimonious yet accurate model we can also employ the Akaike information criterion (AIC) [Akaike1998InformationTheory, Akaike1974StatisticalModelIdentification] and Bayesian information criterion (BIC) [Schwarz1978EstimatingDomensionOfAModel], as in [Mangan2017ModelSelection]. It is also possible to sweep through the parameter and candidate functions , and then choose the Pareto optimal model from a family of models on the Pareto front balancing accuracy and efficiency; this is the approach in the original SINDy work [Brunton2016SINDy] and in earlier work leveraging genetic programming to discover dynamics [Bongard2007pnas, Schmidt2009science]. In this work, we take a different approach, selecting models based on performance on a test data set that has been withheld for model validation. One error function is the model fit:

(11)

In practice, for rational dynamics, we select based upon the predicted derivative :

(12)

For implicit dynamics where each state derivative may be written as a rational function

(13)

then we restrict the candidate functions to for some to identify a separate sparse model for each . Several candidate functions may provide accurate and sparse models. These different models may further be cross-references to check that the same terms are being selected in each model, providing additional information for model selection and validation.

3.2 Constrained Optimization Formulation

Figure 2: Schematic illustrating the constrained formulation of the SINDy-PI algorithm.

In (10) each candidate function was tested individually in a parallel optimization. However, each of these individual equations may be combined into a single constrained system of equations

(14)

We constrain to have zero entries on the diagonal, as shown in Fig. 2, which is the same as removing the candidate function from the library in the separate optimization problems in (10). Without this constraint, the trivial solution will provide the sparsest and the most accurate model. This may be written as a formal constrained optimization problem:

(15)
s.t.

This optimization is non-convex, although there are many relaxations that result in accurate and efficient proxy solutions [Tibshirani1994, Brunton2016SINDy, ZhengSR3]. In this work, we will use sequentially thresholded least squares, so that any entry will be set to zero; the sparsity parameter is a hyperparameter, and each column equation may require a different parameter . The constrained formulation in (15) can be solved efficiently in modern optimization packages, and we use CVX [CVX_Software, Grant_Boyd2008_Nonsmooth_Convex_Programs]. After solving (15) we have numerous candidate models, one for each column of , given by

(16)

The sparse models that result in an accurate fit are candidate implicit models, and they may be assessed using the model selection approaches outlined above. These various models may be cross-referenced for consistency, as the same models will have the same sparsity pattern. This information can then be used to refine the library , for example to only include the nonzero entries in the sparse columns of .

3.3 Noise Robustness

We now compare the noise sensitivity of SINDy-PI and implicit-SINDy on the simple one-dimensional Michaelis–Menten model for enzyme kinetics [Johnson2011MMK_Translate, Mangan2016ImSINDY, Cornish_Bowden2015MMK], given by

(17)

where denotes the concentration of the substrate, denotes the influx of the substrate, denotes the maximum reaction time, and represents the concentration of half-maximal reaction. We use the same parameters as in [Mangan2016ImSINDY], with , , and . Figure 3 shows the result of the noise robustness of SINDy-PI and implicit-SINDy. In this example, SINDy-PI is able to handle over more measurement noise than implicit-SINDy, while still accurately recovering the correct model. Details are provided in App. A.

Figure 3: Performance comparison of SINDy-PI and implicit-SINDy on Michaelis-Menten kinetics (17). The derivative is computed by the total-variation regularization difference (TVRegDiff) [Chartrand2011TVReg] on noisy state measurements. The violin plots show the cross-validated distribution of the number of incorrect terms in the model. The green region indicates no structural difference between the identified model and the ground truth model.

3.4 Data Usage

The data required to correctly identify a model is a critical aspect when comparing SINDy-PI and implicit-SINDy. Many experimental data sets are limited in volume, and thus our goal is to identify a model with as little data as possible. In this section, we compare the SINDy-PI and implicit-SINDy methods on the challenging yeast glycolysis model [Schmidt2011YeastGlycolysis, Mangan2016ImSINDY] given by

(18a)
(18b)
(18c)
(18d)
(18e)
(18f)
(18g)

Equation (18f) is the most challenging equation to discover in this system, and Fig. 4 compares the success rate of SINDy-PI and implicit-SINDy in identifying this equation. SINDy-PI uses about times less data than the implicit-SINDy when identifying (18f). Details are provided in Appendices B and D.

Figure 4: Success rate of SINDy-PI and implicit-SINDy identifying yeast glycolysis (18f) with different lengths of training data.

3.5 Comparison for Implicit PDE Identification

Figure 5: Comparison of SINDy-PI and PDE-FIND on an implicit PDE problem given by the modified KdV equation (19). As we increase , the rational term begins to play a significant role in the system behavior. For small , PDE-FIND compensates for the effect of the rational term by tuning the other coefficients. When is large, PDE-FIND overfits the library. SINDy-PI, on the other hand, correctly identifies the rational term.

We now investigate the ability of SINDy-PI to discover a PDE with rational terms, given by a modified KdV equation

(19)

where is a loss term and is a gain term. We fix and vary the value of from to . As increases, the implicit term gradually dominates the dynamics. Figure 5 shows the results of SINDy-PI and PDE-FIND [Rudy2017PDE-find] for different values of . For large , SINDy-PI is able to accurately identify the rational function term, while this is not possible for PDE-FIND, since this term is not in the library. Details of the identification process are given in App. C.

4 Advanced Examples

We will now demonstrate the SINDy-PI framework on several challenging examples, including the double pendulum, an actuated single pendulum on a cart, the Belousov-Zhabotinsky PDE, and the identification of conserved quantities. All examples are characterized by rational nonlinearities, and we were unable to identify them using SINDy or implicit-SINDy, even in the absence of noise.

4.1 Mounted Double Pendulum

Figure 6: Schematic illustration of SINDy-PI identifying a mounted double pendulum system.

In our first example, we use SINDy-PI to discover the equations of motion of a mounted double pendulum, shown in Fig. 6. The double pendulum is a classic example of chaotic dynamics [Graichen2007], and was an original challenging example used to demonstrate the capabilities of genetic program for model discovery [Schmidt2009science]. Correctly modeling the nonlinear dynamics is vital for accurate control [Graichen2007].

We simulate the double pendulum dynamics, derived from the Euler-Lagrange equations, and use SINDy-PI to re-discover the dynamics from noisy measurements of the trajectory data. The governing equations and SINDy-PI models are provided in App. F. Because these dynamics have rational nonlinearities, the original SINDy algorithm is unable to identify the dynamics, making this a challenging test case. The state vector is given by , and the parameters of the simulation are given in App. D. The training data is generated from an initial condition , simulated for seconds using a time step of seconds. The validation data is generated from an initial condition , simulated for seconds with time step seconds.

To test the robustness of SINDy-PI, we add Gaussian noise to both the training and validation data. We test the resulting models using a new testing initial condition . We construct our library to include over 40 trigonometric and polynomial terms. The most challenging part of this example is building a library with the necessary terms, without it growing too large. The library cannot be too extensive, or else the matrix becomes ill conditioned, making it sensitive to noise. To reduce the library size, we use one piece of expert knowledge: the trigonometric terms should only consist of and , the rotational angles of the pendula.

The candidate functions are chosen as a combination of state derivatives and trigonometric functions. Fig. 6 shows that SINDy-PI can identify the equations of motion for low noise. For larger noise, SINDy-PI misidentifies the dynamics, although it still has short term prediction ability.

4.2 Single Pendulum on a Cart

Figure 7: SINDy-PI is used to identify the single pendulum on a cart system. Control is applied to the cart, and both the cart and pendulum states are measured. When the measurement noise is small, SINDy-PI can identify the correct structure of the model.

We now apply SINDy-PI to identify a fractional ODE problem with control input, given by the single pendulum on a cart in Fig. 7. SINDy has already been extended to include control inputs [Kaiser2018prsa], although the original formulation doesn’t accommodate rational functions.

The dynamics are derived from the Euler-Lagrange equations. All system parameters except for gravity are chosen to be , as summarized in App. D; the governing equations and SINDy-PI models are shown in App. E. The cart position is denoted by . The state vector is given by . The equations of motion are given by

(20a)
(20b)
(20c)
(20d)

Eq. (20) is simulated with a time step of to generate the training and testing data for model selection. The training data is generated using an initial condition with the control input chosen as , for time to . Similarly, the validation data is generated using an initial condition with the control input chosen as , for time to .

The library is constructed using a combination of trigonometric and polynomial terms. Around different basis functions are used for the library, and around terms are tested as candidate functions. We add Gaussian noise to all system states. We then test the SINDy-PI model on a testing initial condition with control input for time to . Fig. 7 shows the resulting SINDy-PI models. The structure of the model is correctly identified up to a noise magnitude of . Beyond this noise level, the SINDy-PI identified model only has short term prediction ability.

4.3 Belousov–Zhabotinsky Reaction

Figure 8: SINDy-PI is able to identify the dynamics of Belousov–Zhabotinsky reaction.

We now apply SINDy-PI to a challenging PDE with rational nonlinearities, the BelousovZhabotinsky (BZ) reaction. The model of the BZ reaction is given by [Vanag2004BZ_Reaction]

(21a)
(21b)
(21c)
(21d)

where , , , and are dimensionless variables and denotes the Laplacian operator.

The strong coupling dynamics and implicit behavior in (21a) make the data-driven discovery of the BZ reaction challenging when using implicit-SINDy and PDE-FIND. However, SINDy-PI correctly identifies the dynamics of BZ-Reaction, as shown in Fig. 8. To generate the data, we use a spectral method [Trefethen_Book, kutz2013data] with time horizon and time step of . We use discretization points with spatial domain ranging from to . The initial condition is chosen to be a mixture of Gaussian functions. of the data is used for training, and the remaining is used for model selection. The right-hand side library is normalized during the sparse regression process. A range of sparsity parameters are tested from to , with increments of The other system parameters in (21) are given in App. D and the SINDy-PI model is given in App. G.

4.4 Extracting Physical Laws and Conserved Quantities

In this final example, we demonstrate how to use SINDy-PI to extract governing physical laws and conserved quantities from data. Many systems of interest are governed by Hamiltonian or Lagrangian dynamics. Instead of identifying the ODE or PDE equations of motion, it might be possible to extract the physical laws directly. These equations contain important information about the system and may be more concise, useful, and straightforward than the underlying ODE or PDE. For example, given a Lagrangian, we can derive the equations of motion.

The most difficult aspect of using SINDy-PI to identify a physical law is how to build the library. Conservation laws may contain higher-order derivatives, such as . To include all possible terms, the library may become exceedingly large. The library size will also increase if the system has many states. Large libraries make the sparse regression sensitive to noise. Thus, extracting the physical law from data using SINDy-PI is still challenging due to the lack of constraints when constructing the library function. We only show one example in our paper to demonstrate that it is possible to achieve this using SINDy-PI, but further work is required to reduce the library size so that the sparse regression is robust.

Figure 9: SINDy-PI is used to extract the conserved quantity for a double pendulum.

As an example, we consider the double pendulum shown in Fig. 9, with the system parameters given in App. D. In this case, we also account for the friction in the pendulum joint, with friction constants of and for the pendulum arms, respectively. In this case, we extract the Lagrangian of the double pendulum [Graichen2007] using SINDy-PI. To extract this Lagrangian, we simulate the system with initial condition from to with time step . The resulting model is shown in Fig. 9.

5 Conclusions and Future Work

In this paper, we develop SINDy-PI (parallel,implicit), a robust variant of the SINDy algorithm to identify implicit dynamics and rational nonlinearities. SINDy-PI overcomes the sensitivity of the previous implicit-SINDy approach, which is based on a null-space calculation, making it highly sensitive to noise. Instead, we introduce both parallel and constrained optimizations to test candidate terms in the dynamics, making the new SINDy-PI algorithm as robust as the original SINDy algorithm. We also extend the algorithm to incorporate external forcing and actuation, making it more applicable to real-world systems. We demonstrate this approach on several challenging systems with implicit and rational dynamics, including ODEs, actuated systems, and PDEs. In particular, we discover the implicit dynamics for the BZ chemical reaction PDE, the double pendulum mechanical system, and the yeast glycolisis model, which have all been challenging test cases for advanced identification techniques. Throughout these examples, we demonstrate considerable noise robustness and reductions to the data required, over the previous implicit-SINDy algorithm.

Despite the advances outlined here, there are still many important avenues of future work. One limitation of this approach, and of SINDy in general, is in the design of the library of candidate functions. The goal is a descriptive library, but the library size grows rapidly, which in turn makes the sparse regression ill-conditioned. Recently, tensor approaches have been introduced to alleviate this issue, making libraries both descriptive and tractable [Gelss2019mindy], and this is a promising approach that may be incorporated in SINDy-PI as well. More generally, automatic library generation, guided by expert knowledge, is an important topic. Other research directions will involve parameterizing elements of the library, so that the algorithm simultaneously identifies the model structure and the parameters of the sparsely selected terms. Recent unified optimization frameworks, such as SR3 [ZhengSR3, Champion2019SINDySR3], may make this possible. Model selection is another key area that will required focused attention. Balancing accuracy on test data, sparsity of the model, and the potential for overfitting are all serious concerns. The sparse regression and optimization may also be improved for better noise robustness. Finally, modifying SINDy-PI to incorporate prior physical knowledge and to only model the discrepancy with an existing model [KK2019Discrepancy] will be the focus of ongoing work.

Acknowledgments

SLB acknowledges support from the Army Research Office (ARO W911NF-19-1-0045) and the Air Force Office of Scientific Research (AFOSR FA9550-18-1-0200). JNK acknowledges support from the Air Force Office of Scientific Research (AFOSR FA9550-17-1-0329). We also acknowledge valuable discussions with Aditya Nair, Eurika Kaiser, Brian DeSilva, Tony Piaskowy, Jared Callaham, and Benjamin Herrmann. We also thank Ariana Mendible for reviewing the manuscript and providing useful suggestions.

References

Appendix A Noise Sensitivity of SINDy-PI and implicit-SINDy

a.1 Performance Evaluation Criteria

To compare the performance of SINDy-PI and implicit-SINDy for noisy data, we must define an evaluation criteria. We compare the performance of the best model generated by each method that has the lowest prediction error on the test data, selected according to Eq. (12). To compare the models generated by the two methods with the ground truth model, we use the concept of model discrepancy [KK2019Discrepancy, M.C.Kennedy2001BayesianCalibration, Arendt2012ModelUncertainty] and set prediction accuracy, structural accuracy, and parameter accuracy as our performance criteria. A good prediction error does not guarantee the model has good structural accuracy and parameter accuracy, and vice versa, motivating multiple performance criteria.

a.2 Numerical Experiments

We use the Michaelis–Menten kinetics, given by Eq. (17), to compare the performance of SINDy-PI and implicit-SINDy. We performed our numerical experiments as follows:

  1. Randomly generate 2400 different initial conditions of different magnitudes ranging from to . Simulate those initial conditions using a fourth-order Runge-Kutta method with time step and time horizon .

  2. Add Gaussian noise to the data. different Gaussian noise instances with magnitudes ranging from to are used.

  3. Compute the derivative of the noisy data. We investigate several approaches, including finite-difference and total-variation regularized difference (TVRegDiff) [Chartrand2011TVReg] derivatives. In all cases, SINDy-PI is several orders of magnitude more robust to noise than implicit-SINDy. TVRegDiff generates more accurate derivatives, but also requires hyperparameter tuning and causes aliasing, so we trim the ends of the time series (first and last ). It is possible to add Gaussian noise to the clean derivative data to investigate robustness, although this is less relevant for real-world scenarios, where only noisy state data is available.

  4. Train SINDy-PI and implicit-SINDy models on noisy training data. For each noise level, we sweep through different sparsity parameters for SINDy-PI, from to . The is varied by a factor of 2 [Qu2014ADM] to calculate the null space in the implicit-SINDy method. The library for implicit-SINDy and SINDy-PI is

    (22)
  5. Due to the various parameter values, we use model selection to choose a model. We use the test data with the same noise magnitude to perform the model selection process. The ratio of training data and testing data is .

  6. The best model generated by the two methods are compared. We use the prediction error, error in the model structure, and parameter error as our model performance evaluation criteria. We average the performance over random initial conditions.

Many parameters affect the performance of these methods: the length of training data, prediction steps to calculate prediction error, the initial conditions for training data, choice of the library, and the derivative computation. We have attempted to carefully optimize each method, although an exhaustive parameter sweep is beyond the scope of the present work. However, in all cases SINDy-PI outperforms implicit-SINDy.

Appendix B Data Usage of SINDy-PI and implicit-SINDy

Sec. 3.4 investigates the data usage of SINDy-PI and implicit-SINDy on the yeast glycolysis model in Eq. (18). The parameters of this problems are given in Table. 1. The data usage comparison is performed by the following steps:

  1. Generate training data by simulating Eq. (18) with parameters in Table. 1 and a time step of , with time horizon . We simulate the system using random initial conditions with magnitude ranging from to .

  2. Randomly shuffle the training data and select the first percent to train SINDy-PI and implicit-SINDy models. No noise is added since we only care about the effect of the data length in this case. The sparsity parameter is fixed for both algorithms (different values); this value is selected for a single percentage where both methods fail to identify the correct model, and we sweep through .

  3. Run the numerical experiment for times for each data length and calculate the percentage of times the two algorithms yield the correct structure of the Eq. (18f).

Parameter
Value 2.5 -100 13.6769 200 13.6769 -6 -6 6 -64 6 16 64 -13
Parameter
Value 13 -16 -100 1.3 -3.1 -200 13.6769 128 -1.28 -32 6 -18 -100
Table 1: The parameter used for simulating the Eq. (18).

The final comparison is shown in Fig. 4. Data usage requirements for other state equations are given in Table  2; Fig. 4 shows results for the hardest equation to identify. The other equations require less data. Normalizing the SINDy-PI library improves data learning rates as well.

Equation Eq. (18a) Eq. (18b) Eq. (18c) Eq. (18d) Eq. (18e) Eq. (18f) Eq. (18g)
Library Order
SINDy-PI un-normalized Left-Hand Side
Threshold
Data Usage
SINDy-PI normalized Left-Hand Side
Threshold
Data Usage
implicit-SINDy normalized Threshold
Data Usage
Table 2: Comparison of data usage of SINDy-PI and implicit-SINDy on other states.

Appendix C SINDy-PI and PDE-FIND on Rational PDE Problem

In Sec. 3.5, we compared the performance of SINDy-PI and PDE-FIND on a modified KdV equation. The simulation data is obtained using a spectral method [Trefethen_Book] with a time step of and time horizon , spatial domain , and spatial discretization points. The library of PDE-FIND is chosen to be

(23)

and the right-hand side library for the SINDy-PI is chosen to be

(24)

while the left-hand side library is chosen to be

(25)

For both SINDy-PI and PDE-FIND, we used different values for the sparsity parameter ranging from to with step size . We use of the simulation data for training and for testing and model selection. We calculate the normalized prediction error for all models on state and the model with minimum prediction error is selected as the final model.

Appendix D Parameter Values for Simulations

The parameters for the double pendulum simulation in Sec. 4.1 are given in Table. 4. The parameters used to simulate the Belousov-Zhabotinsky reaction model in Eq. (21) are given in Table. 5.

Parameter
Value
Table 3: Parameters used to simulate the double pendulum.
Parameter
Value
Table 4: Parameters used to simulate the single pendulum on a cart.
Parameter
Value
Table 5: Parameters Used in Eq. (21) for Simulating the Belousov-Zhabotinsky Reaction Model.

Appendix E SINDy-PI Models for the Single Pendulum on a Cart

The Lagrangian for the single pendulum on a cart with an input force on the cart is:

(26)

where is the mass at the end of the pendulum arm, is the mass of the cart, is the length of the pendulum arm, is the position of the cart, and is the pendulum angle. We do not consider damping in this case. Using the numeric values and this simplifies to

(27)

The Euler-Lagrange equation of the system are

(28a)

where is the force applied to the pendulum cart. It is possible to isolate and :

(29a)
(29b)

It is possible to write this as a system of four coupled first-order equations

(30a)
(30b)
(30c)
(30d)

With the numerical values used here, this becomes

(31a)
(31b)
(31c)
(31d)

With no noise added, SINDy-PI discovers the correct equations:

which is correct.

When we add random noise with magnitude of SINDy-PI discovers

which is the correct structure with small parameter errors.

When we add random noise with magnitude of the SINDy-PI incorrectly identifies

When we add random noise with magnitude of the SINDy-PI incorrectly identifies

Appendix F SINDy-PI Models for the Mounted Double Pendulum

For a mounted double pendulum system shown in Fig. 6 we could have following parameters: the parameters of the pendulum are center of mass and , center of mass position and , arm length and , arm inertia and , arm rotational angle and , gravity acceleration . Those values could be seen from Table. 4. If we consider friction between the pendulum joint, we could define and as our damping coefficient.It is easy to derive the Lagrangian of the mounted double pendulum which is given by

(32)

The damping term caused by friction with friction coefficients and is

(33)

The Euler-Lagrange equations with a Rayleigh dissipation term are then:

(34a)
(34b)

The symbolic form of the Eq. (34a) is

(35)

and the symbolic form of Eq. (34b) is