Sparse identification of Lagrangian for nonlinear dynamical systems via proximal gradient method

The autonomous distillation of physical laws only from data is of great interest in many scientific fields. Data-driven modeling frameworks that adopt sparse regression techniques, such as sparse identification of nonlinear dynamics (SINDy) and its modifications, are developed to resolve difficulties in extracting underlying dynamics from experimental data. However, SINDy faces certain difficulties when the dynamics contain rational functions. The Lagrangian is substantially more concise than the actual equations of motion, especially for complex systems, and it does not usually contain rational functions for mechanical systems. Few proposed methods proposed to date, such as Lagrangian-SINDy we have proposed recently, can extract the true form of the Lagrangian of dynamical systems from data; however, these methods are easily affected by noise as a fact. In this study, we developed an extended version of Lagrangian-SINDy (xL-SINDy) to obtain the Lagrangian of dynamical systems from noisy measurement data. We incorporated the concept of SINDy and used the proximal gradient method to obtain sparse Lagrangian expressions. Further, we demonstrated the effectiveness of xL-SINDy against different noise levels using four mechanical systems. In addition, we compared its performance with SINDy-PI (parallel, implicit) which is a latest robust variant of SINDy that can handle implicit dynamics and rational nonlinearities. The experimental results reveal that xL-SINDy is much more robust than the existing methods for extracting the governing equations of nonlinear mechanical systems from data with noise. We believe this contribution is significant toward noise-tolerant computational method for explicit dynamics law extraction from data.


I. INTRODUCTION AND RELATED WORKS
Since the early modern history of humanity, scientists have always been trying to come up with models that can capture real-world phenomena.Such models are desired because they can be used to devise solutions to real-world problems.For centuries, the process of refining hypothesis and models from observation data have been conducted manually.Automating this process has long been of great interest in the scientific community.
Many attempts have been made to extract physical laws autonomously from data.With the abundance of data and cheaper yet powerful hardware, the deep learning-based methods have gained a lot of attraction and have been widely used to model and control dynamical systems [1], [2].It has also been shown that deep learning is also capable of approximating invariant quantities from dynamical systems such as the Hamiltonian [3] and the Lagrangian [4].However, deep learning models act as black boxes; it does not provide The authors are with the Neuro-Robotics Lab, Department of Robotics, Tohoku University, 980-8579, Sendai, Japan (email : adam.syammas.zaki.purnomo.p8@@dc.tohoku.ac.jp; hayashibe@tohoku.ac.jp) insights into how each observation variable affects and relates to each other.
Recent trends, however, favor parsimonious models, models with the lowest complexity to describe the observation data.A ground-breaking work done by Schmidt and Lipson [5] shows us that it is possible to extract the governing mathematical expressions from observation data.Symbolic regression is used to find the nonlinear differential equations that describe the behavior of the system, but symbolic regression tends to be expensive.The sparse identification of nonlinear dynamics (SINDy) [6] models nonlinear differential equations of dynamics as a linear combination of nonlinear candidate functions and obtains a parsimonious model through sparse regression.
While there are many applications of SINDy across different fields [7]- [10], SINDy faces certain difficulties when the dynamics contain rational functions.Including rational functions into the library of candidate functions would tremendously increase the size of the library, making the sparse regression challenging.A modification of SINDy, implicit-SINDy [11], reformulates the SINDy problem into implicit form to address this challenge, albeit this method is sensitive to noise.SINDy-PI [12] is proposed to improve the performance of implicit-SINDy in terms of noise robustness.While SINDy-PI is much more robust than implicit-SINDy, it can only obtain the correct dynamical structure with noise magnitude on the scale of up to 10 −3 which might not be sufficient for a real-world application.On top of that, if the incorrect combination of denominator terms is discovered, the predicted system might blow up when the denominator is equal to zero The principle of least action is fundamental to many dynamical systems [13].The principle states the trajectory chosen by the system is the one that minimizes a certain cost function.This cost function is the so-called 'action,' which is defined as the integral of the Lagrangian for an input evolution between a certain time period.Compared to the underlying differential equations, Lagrangian has a desirable property in which it is a single scalar quantity that contains all information to predict the behavior of the systems.In robotics, the derivation of the dynamics often starts from the Lagrangian of the systems.
Several works have been proposed to approximate the Lagrangian from data with polynomial basis functions [14], [15].However, approximating Lagrangian with polynomial basis functions will only be useful for a particular trajectory of the system and will not likely generalize well across different initial conditions.Lagrangian-SINDy [16] is a SINDy-based method designed to extract the Lagrangian of nonlinear dynamics and is shown to be able to retrieve the true form of Lagrangian of several dynamical systems.However, the author mentioned in the paper that Lagrangian-SINDy is sensitive to noise and cannot recover the Lagrangian when the training data is corrupted by Gaussian noise even with magnitude in the scale of 10 −7 .Noise will always present in real-world systems, and developing a robust method against noise is important for real-world applications.
In this work, we propose a method called extended Lagrangian-SINDy (xL-SINDy) that can discover the true form of the Lagrangian and is more robust in the presence of noise compared to Lagrangian-SINDy and SINDy-PI.We demonstrated the effectiveness of xL-SINDy against different noise levels and compare the robustness of xL-SINDy against SINDy-PI in physical simulations with four dynamical systems: A single pendulum, a cart-pendulum, a double pendulum, and a spherical pendulum.This paper is organized as follows; section II describes how we use basic ideas from previous works to formulate the problem and develop the learning method, section III presents the results of simulation experiments on the aforementioned dynamical systems, and section IV provides closing and remarks.

A. Problem Formulation
Inspired by the concept of SINDy [6], we consider a Lagrangian expression in a structure of a linear combination of nonlinear candidate functions.Let q = (q 1 , q 2 , ..., q n ) be the configuration of a system in a generalized coordinate of a system, the Lagrangian of the system is expressed as where, φ k (q, q), k = 1, ..., p are a set of nonlinear candidate functions, and c k , k = 1, ..., p are the corresponding coefficients.We are interested to find the value of c = (c 1 , c 2 , ..., c p ) where we believe that the majority of the coefficients are zero.The Lagrangian of the system satisfies the Euler-Lagrange equations given by where (∇ q ) i ≡ ∂ ∂qi .We consider three different scenarios: • Case I: External input τ ext of the system is provided.
• Case II: No external input is provided.
• Case III: Prior Lagrangian knowledge of a simpler system that forms a constituent of the system is provided.1) With External Input: In the case where input τ ext is provided, substituting (1) in (2) yields where τ pred is the predicted value of the external input τ ext given a set of coefficient c = (c 1 , c 2 , ..., c p ).We can further expand the time derivative d dt by using chain rule, giving us the terms q and q To avoid verbose notation, we define the following notations Substituting ( 5), (6), and ( 7) in (4) yields and we define the following cost function that we want to minimize to obtain the Lagrangian of the system We can use the above cost functions if the external torque is provided.In the case of passive systems where no external torque is provided, we will just end up minimizing the residual cost function J(c) = −τ pred (c) 2  2 .From equation (8), we can observe that the τ pred (c) is in the form of linear combination of coefficient c.Thus, minimizing this residual cost function is equivalent to the problem of finding a sparse null space which is an arduous task with current optimization methods.This will lead us to the formulation of a new cost function that will be explained in the next section.
2) Without External Input: In the case of passive systems, in which no external input τ ext is not provided, eq. ( 8) can be modified so that we can solve for qpred expressed as where qpred represents the predicted value of acceleration q and (•) −1 represents matrix inverse.In practice, we use Moore-Penrose pseudo inverse to calculate eq.(10) to avoid numerical instability.We define the following cost function to learn the Lagrangian of the system Due to the inverse operation, equation ( 11) is non-convex with respect to variable c, making the optimization process not always converge to the global minimum.We empirically found that with a library of more than 20 candidate functions, the learning process will hardly converge, and the value of the cost function rarely touches a single-digit value even after a long period of iterations.Therefore, we only use this case only when no prior knowledge is available and the system is not complex such as a single pendulum.For more complex systems, such as a multi-DOF system, it is preferable that the external input τ ext is provided.In the case where no external input is provided, prior Lagrangian knowledge of a simpler system that forms a constituent of the larger system can be utilized to boost the learning process which will be explained in the next section.
3) With Prior Knowledge: For multi-DOF non-relativistic systems, the Lagrangian can be described as , where T i and V i are the kinetic energy and potential energy of each constituent of the system.Since the total Lagrangian of the system is the sum of the Lagrangian of its constituents, it is reasonable to assume that the nonlinear terms that appear in each constituent will also appear in the total Lagrangian of the system [16].
Given prior knowledge of a constituent of the system, we pick one out of several terms that appear in the total Lagrangian of the system and label them as φ r (q, q).The Lagrangian of a system is not unique; many forms of Lagrangians can satisfy Euler-Lagrange's equation for a particular system.For example, L ′ = kL, where k is a constant, still satisfies the Euler-Lagrange's equation.By multiplying equation (1) with k = 1 cr , eq. ( 1) can be modified as where c ′ k = c k cr .Now, the variable c ′ k becomes the coefficients that we are interested to solve.We can simply redefine c k := c ′ k for notation simplicity.The Euler-Lagrange equation of the system can be expressed as We define the following notation where Υ lef t and Υ right represent the left hand side (LHS) and the right hand side (RHS) of eq. ( 13).By minimizing the following cost function it is possible to obtain the true Lagrangian of the system.We usually have more than one option of φ r to construct Υ right .In practice, we have to test all of them one by one and choose the one that yields the best model.

B. Learning Lagrangian
The proposed learning method to obtain the Lagrangian is summarized in Fig. 1.We start with the problem formulation as described in the previous section.Given a dynamical system, we gather times series data {t i , q(t i ), q(t i ), q)(t i , τ ext (t i )} N i=1 from several initial conditions.We then proceed to construct a library of candidate functions.
In general, the larger the library of the candidate functions, the more difficult the optimization problem becomes.It is especially true when several candidate functions can behave in a similar manner, such as in the trigonometric family functions.It is important to carefully construct a sufficient library but not too large so that the optimization problem is still tractable.It is also better, however, not to include trivial terms that satisfy Euler-Lagrange's equation regardless of trajectories such as L = q n q.Technically, even though trivial terms will not affect the behavior of the systems, it is better to remove them so that we can reduce unnecessary complexity in our library.Depending on the case of the problem, a different cost function should be defined.
As mentioned previously, we believe that the correct solution is sparse where the majority of the coefficients are zero.Therefore, we add the L1 regularization term to the cost function for sparsity constraint [17] expressed as where λ is the sparsity promoting parameter that we have to carefully tune.In this work, we use the accelerated proximal gradient descent method [18] to minimize the composite cost function defined above.Given an initial point c 0 , the update step of proximal gradient descent is defined as where c i is the coefficient c at iteration i, α is the learning rate, and prox L1 (•) is the proximal operator for L1 norm.The L1 norm penalty term is a separable sum of the component of its input, and a proximal operator is used to minimize this term.The proximal operator for the L1 norm is well defined separately for each component of the input and expressed as follows, where k is the k th entry of the input vector β.As for case III, if we know other terms that appear in the Lagrangian but are not used to construct Υ right , we don't put a penalty on these terms by not applying the proximal operator in eq. ( 20) for index k corresponding with these terms.We proceed to initialize the value of the coefficient c, the learning rate α, and the L1 norm penalty parameter λ.The learning process is done in several stages, with 100 epochs and batch size equal to 128 for each stage, until the cost function reaches the defined tolerance value as shown in Fig. 1.The tolerance value is ideally at 10 −3 .However, converging to this value might not be possible in the presence of noise, and we have to relax the tolerance value; otherwise, the algorithm will never stop.In the beginning, the number of candidate functions in the library is usually large, and we want to eliminate non-relevant candidate functions as much as possible during the first learning stage.Therefore, we initially set the value of λ to be quite high, which is between 1 and 5.
It is important to note that every candidate function may have different magnitude scales.The L1 norm penalty penalizes all terms equally regardless of the magnitude scale, resulting in candidate functions with smaller magnitude scales being penalized more.It may or may not be necessary to do scaling in the first learning stage, where the value of λ is high, by multiplying each candidate function with scaling term s k in eq. ( 1) for cases I and II, or eq.( 12) for case III depending on the differences of magnitude scale between each candidate function.The learning rate is also an important hyper-parameter, especially during the first learning stage.We found that a high learning rate in the Fig. 2: Dynamical systems used to verify xL-SINDy.From upper left to bottom right: A single pendulum, a cart pendulum, a spherical pendulum, and a double pendulum.For all systems, the length of the rod is L = 1.0 m, the mass of all pendulums are m = m p = m 1 = m 2 = 1.0 kg, and the gravitational acceleration is g = 9.81 m/s 2 .For the cart-pendulum, the mass of the cart is m c = 0.5 kg.initial stage can cause the relevant terms to be penalized, preventing the model from obtaining the true Lagrangian of the system.During the initial stage, we set the learning rate α <= 10 −5 .
At the end of every learning stage, we perform hardthresholding by removing index k from eq. ( 1) or (12), where the value of c k < threshold.This step effectively reduces the number of candidate functions considered in the learning process, making the convergence much faster than if hardthresholding were not performed.We then check whether the cost function has reached the tolerance or not.If it is the latter, we proceed to the next learning stage.With fewer candidate functions after the previous hard-thresholding process, we can decrease the value of λ and increment the learning rate α to speed up the learning process.The rate of how we increment α and decrement λ also matters.If we increment α too fast, the optimization process can even take longer to converge, and if we decrement λ too fast, we might not be able to get rid of all of the non-relevant terms.A tuning process is needed to find the most appropriate rate to increase α and decrese λ.This step is repeated over and over again until the cost function reaches the tolerance value.Normally, the learning process will take around 3 or 4 stages until the tolerance value is reached.Once the tolerance value is reached, we compute the value of the coefficient to the eq.( 1) for cases I and II, or (12) for case III, and we obtain the analytical form of Lagrangian of the system.

C. Dynamical Systems and Experiments
We evaluated xL-SINDy with four ideal dynamical systems as shown in Fig. 2. In this work, at first, we test xL-SINDy when the dynamical systems are excited by external inputs τ ext = f (sin ωt, cos ωt), where ω is a random frequency.The computation used to learn the model is described by the computation of case I.
We also test the case of passive systems where no external input τ ext is provided.For the case of a single pendulum, we assume that no prior knowledge is available.Thus, we use the computation described in case II.For the cart pendulum, double pendulum, and spherical pendulum, these systems either have a single pendulum as one of their constituents or a more complex version of a single pendulum.Thus, assuming that we already obtained the Lagrangian of a single pendulum, we can use the computation described in case III to bootstrap the learning process.
For each system, we collect training data by performing simulation with 100 initial conditions for a period of 5s each and 100 Hz of measurement frequency.After obtaining the analytical form of the Lagrangian, we create a validation data set to test the obtained model by calculating the predicted states for the accuracy evaluation.We compute the Euler-Lagrange's equation with the obtained model, retrieve the differential equation of the system, and integrate the equations to compare them with the actual validation data.We also tested our method with training data that are corrupted by zero-mean white Gaussian noise N (0, σ) on different scale magnitude in the range of 10 −8 <= σ <= 10 −1 .
Finally, we also compare the performance of xL-SINDy on several passive dynamical systems with noisy training data against SINDy-PI [12].

III. RESULTS
In this section, we demonstrate the effectiveness of xL-SINDy with the aforementioned nonlinear dynamical systems against various noise levels.The obtained Lagrangian for each system is summarized in Table.I for the case of active systems, and in Table.II for the case of passive systems.The performance of xL-SINDy against the true model from our simulation experiments on a cart pendulum, a double pendulum, and a spherical pendulum is shown in Fig. 3 in the case of active systems, and in Fig. 4 in the case of passive systems along with the comparison against SINDy-PI.From Table.I and Table.II, it can be seen that xL-SINDy performs better to extract the correct structure if external inputs are provided.It includes fewer wrong additional terms in the model compared to when no external input is provided.
As for the performance of xL-SINDy against SINDy-PI, in the second column of the plot in Fig. 4, where the noise magnitude is σ = 2 × 10 −2 , we can observe that SINDy-PI has already started to deviate from the true models in all three dynamical systems.At the same noise magnitude, xL-SINDy still predicts accurate models.It is also good to note that the model estimate of xL-SINDy is still reasonable even though wrong additional terms are included in the Lagrangian from the example of the cart pendulum under the noise magnitude of σ = 6 × 10 −2 .It indicates that the model estimate is potentially usable even when an incorrect Lagrangian structure is discovered.Our simulation results demonstrate that xL-SINDy has notably better prediction accuracy compared to SINDy-PI in the presence of higher noise magnitude in all three dynamical systems.

A. Single Pendulum
The state of a single pendulum is described by [θ, θ], and the Lagrangian expression of a single pendulum is given by L = 1 2 m θ2 + mg cos θ, Substituting the parameter given in Fig. 2, the true Lagrangian expression is shown in the second row and second column of Table.I. To construct a library of candidate functions, we create a polynomial combination of {θ, θ, cos θ, sin θ} up to the second order while excluding trivial terms such as θ and θ θ resulting in 12 candidate functions.Training data with initial conditions of [−π < θ < π, 0] are created.
The initial value of the hyperparameters are α = 10 −5 and λ = 0.1.The cut-off threshold is 10 −2 for the initial learning stage and 10 −1 for the subsequent learning stages.In the subsequent learning stages, α is increased by a factor of 2, and λ is decreased by a factor of 10.The training converged in three stages for noise magnitude σ <= 10 −3 , and four stages for higher magnitude with a relaxed tolerance value.
In the case of an active system, the correct Lagrangian structure, the ones without additional terms or missing terms compared to the true Lagrangian form, can be obtained in the presence of noise magnitude up to σ = 1 × 10 −1 .While in the case of a passive system, the correct Lagrangian structure can be obtained with noise magnitude up to σ = 6 × 10 −2 .Even though the coefficients obtained differ from the true model, the ratio of coefficients between the two terms is close compared to the true model.

B. Cart Pendulum
The state of the cart pendulum is represented as [θ, θ, x, ẋ], and the Lagrangian with numerical coefficients is shown in the second row and third column of Table .I with parameters given by Fig. 2. A library of candidate function with a polynomial combination of { θ, cos θ, sin θ, x, ẋ} up to the third order is constructed, resulting in 55 candidate functions.We here exclude the term θ because this term does not appear in the Lagrangian of a single pendulum system.Training data with initial conditions of [−π < θ < π, 0, 0, 0] are created.
The Lagrangian of a single pendulum contains θ2 and cos θ.Hence, both terms will also appear in the Lagrangian of the cart pendulum.We tested both θ2 and cos θ to construct Υ right as described in eq. ( 15), and we found that the term θ2 gives better results.The initial value of the hyperparameters are α = 10 −5 and λ = 1.The cut-off threshold, the increment of α, and the decrement of λ are the same as in the previous case.The training converged in three stages for noise magnitude σ <= 2 × 10 −2 , and four stages for higher magnitude with a relaxed tolerance value.
From the table, xL-SINDy can recover the correct structure of the Lagrangian with noise magnitude up to σ = 1 × 10 −1 with external input and noise magnitude up to σ = 4 × 10 −2 without external input.In contrast, SINDy-PI can only recover the correct structure with noise magnitude up to σ = 5 × 10 −3 .Therefore, in the case of the cart pendulum, xL-SINDy is 8 times more robust than SINDy-PI in the presence of noise.On top of that, with a large magnitude of noises, SINDy-PI sometimes predicts a model which blows up as it is shown in the second column of Fig. 4.This happens because it incorrectly predicts the denominator terms.Unlike SINDy-PI, xL-SINDy deals with the Lagrangian instead of the actual equation of motion.Since there is no denominator in the Lagrangian, xl-SINDy still gives reasonable predictions even though wrong additional terms are included in TABLE II: Extracted Lagrangian from simulation data with various noise levels when no external input is provided.For the single pendulum, the result is obtained with a computation described by case II, while the rest are obtained with a computation described by case III with the knowledge of the Lagrangian of a single pendulum.the Lagrangian.

C. Double Pendulum
Given the state of a double pendulum, [θ 1 , θ 2 , θ1 , θ2 ] and the system parameters in Fig. 2, the expression of the Lagrangian with numerical coefficients is shown in the second row and fourth column of Table.I. To build a library of candidate functions, we first separate the set of trigonometric terms {cos θ 1 , sin θ 1 , sin θ 1 , sin θ 2 }, and the non trigonometric terms { θ1 , θ2 }.For each set, we create a polynomial combination up to the second order, resulting in 14 candidate functions and 5 candidate functions respectively.We then generate cross terms between the two sets creating 70 candidate functions, and we have in total 89 candidate functions in the library.Training data are created with initial conditions of Both constituents of the double pendulum are single pendulums.Hence, we have 4 options to construct Υ right : θ1 used to construct Υ right .The initial value of the hyperparameters are α = 5×10 −6 and λ = 1.The cut-off threshold, the increment of α, and the decrement of λ are the same as in previous cases.From our experiments, xL-SINDy can identify the correct structure with noise magnitude up to σ = 10 −1 in both active and passive cases, while SINDy-PI can extract the correct structure of equations of motions with noise magnitude up to σ = 10 −2 .Hence, xL-SINDy is 10 times more robust against noise than SINDy-PI in this experiment.
As it can be seen from the summary table that xL-SINDy is quite robust in the case of double pendulum compared to other dynamical systems.One possible reason why xL-SINDy is more robust in the case of the double pendulum is due to the chaotic signal caused by the double pendulum.For every initial condition, the double pendulum will yield

D. Spherical Pendulum
The state of the cart pendulum is represented as [θ, φ, θ, φ], and the true Lagrangian expression is displayed in the second row and fifth column of Table .I. As in the case of the double pendulum, we first separate the trigonometric terms {cos θ, sin θ} and the non-trigonometric terms { θ, φ, φ}, create polynomial combinations for both sets up to the second order and add cross terms between the two sets.In total, we have 59 candidate functions in our library.The training data are created with initial conditions of [π/3 < θ < π/2, 0, 0, π].We deliberately choose high value of θ and φ as the initial conditions because the equation of motion contains 1 sin θ which could blow up for small value of θ.The spherical pendulum is a higher dimensional analog of a single pendulum in the first case.Therefore, we can think of the Lagrangian of a spherical pendulum as the sum of the Lagrangian of the pendulum in θ direction and φ direction.Since we already know the Lagrangian of a single pendulum in θ direction, we can use θ2 and cos θ to construct Υ right .The initial value of the hyperparameters are α = 1 × 10 −5 and λ = 1.The cut-off threshold, the increment of α, and the decrement of λ are the same as in all previous cases.In this experiment, xL-SINDy is robust only up to σ = 2×10 −2 for both active and passive system cases.In contrast, SINDy-PI is only robust up to σ = 1x10 −3 .Thus, xL-SINDy shows 20 times more robustness against noise than SINDy-PI.
From the second column of Fig. 3 and Fig. II, it can be seen that even though the correct structure can be obtained for the spherical pendulum at σ = 2×10 −2 , the performance of xL-SINDy in the case of an active system is worse.This is because in the case of the passive system, a Lagrangian multiplied by a constant is still a valid Lagrangian, thus as long as the ratio between each coefficient is the same as the true model, then the obtained model is also correct.However, in the case of the active system, since there is an external constrain, the Lagrangian is unique.So, even though the correct structure is obtained, if the coefficient does not closely match the true model, it will less accurate long-term prediction ability.

IV. DISCUSSION
From the results of simulations, it can be concluded that xL-SINDy is more robust to obtaining the correct Lagrangian structure if external input is provided to the system.However, with larger noise, even though xL-SINDy can obtain the  correct Lagrangian structure, it does not guarantee to provide an accurate long-term prediction ability.This is because in the case where external input is provided, the Lagrangian of the system is unique, thus a mismatch coefficient will result in deviation in the long-term prediction.However, this will not be a problem if there is no external input as the Lagrangian is not unique.As long as the ratio between each coefficient is close to the true model, it can still give a good long-term prediction model.However, with the absence of external input, xL-SINDy is a little bit less robust to find the correct Lagrangian structure.
To compare xL-SINDy with other methods, we use passive cases as the baseline of xL-SINDy.The comparison of xL-SINDy, SINDy-PI, and Lagrangian-SINDy is summarized in Fig. 5.In all three dynamical systems used as a comparison, xL-SINDy outperforms other methods in terms of noise robustness.Our experiment results demonstrate that xL-SINDy can overcome the challenge faced by Lagrangian-SINDy.xL-SINDy is capable of discovering the correct Lagrangian expression for idealized nonlinear dynamical systems in the presence of much higher noise magnitude.On top of that, xL-SINDy successfully extracts the Lagrangian in cases where Lagrangian-SINDy fails to do so, such as the non-actuated spherical pendulum [16].The obtained coefficients may not be precisely the same as the true models, but the ratio between coefficients of each term is close to the true models.
From our experiments, while SINDy-PI is also robust against noise up to a certain magnitude, xL-SINDy is 8 to 20 times more robust against noise.SINDy-PI attempts to seek the expression of the dynamics which may contain rational functions.To do so, SINDy-PI reformulates the problem into implicit form, and it requires the library to include candidate functions of the states and the time derivative of the states of the systems.As a consequence, SINDy-PI also has to include q, q, and q variables to build terms in the library.Unlike SINDy-PI which deals with the actual equation of motion, xL-SINDy deals with Lagrangian that only requires q and q variables to build the terms in the library.Hence, naturally, to obtain the same order of a family function (e.g.second order of polynomial functions), xL-SINDy would contain fewer terms in the library.
The second advantage of xL-SINDy over SINDy-PI is that xL-SINDy can have more coverage of various terms with fewer coefficient parameters in the actual equation of motion.This is due to the inherent nature of the Lagrangian formulation itself.For example, let's consider a Lagrangian that only contains one term given by If we substitute this into Euler-Lagrange's formula, we will get the equation of motion that contains three different terms expressed as As it can be seen from the equation above, even though we have three different terms in the equation of motion, all of them correspond to the same coefficient c 0 .This is not the case with SINDy-PI as it would require three different coefficients to represent different terms in the equations of motion.Hence, with fewer parameters, it would be easier to learn the model with xL-SINDy than SINDy-PI while maintaining the same amount of coverage of possible terms in the actual equations of motion.This is one of the reasons why xL-SINDy is more robust than SINDy-PI.
Another reason why xL-SINDy is more robust than SINDy-PI is how the learning process is done.SINDy-PI uses the sequential threshold least-square method [12] where the basic idea is to run the least square method and remove terms with low coefficient sequentially (hard-thresholding).However, we experimentally found that the sequential leastsquare method often fails to remove non-relevant terms if the training data is corrupted.Instead, we combine the idea of hard-thresholding from the sequential least-square method and soft-thresholding from Lasso regression with the proximal gradient method.We found that it performs better to remove non-relevant terms with a higher level of noise.As it can be seen in the table III, in the case of the double pendulum and spherical pendulum, even though xL-SINDy has more terms in the library than SINDy-PI, xL-SINDy still performs better than SINDy-PI with a higher level of noise.
Finally, SINDy-PI may have a problem when the incorrect combination of denominator terms is discovered.In rational function, when the denominator is equal to zero, its value blows up.Indeed, this is the case of SINDy-PI in the case of the cart pendulum and the spherical pendulum from our experimental results.On the other hand, xL-SINDy only contains rational terms since it is a Lagrangian mechanic function.So, it would minimize the possibility of a model that blows up due to incorrect terms.
Like other learning-based methods, xL-SINDy introduces several hyperparameters during the learning process, such as the sparsity constrain λ, the learning rate α, the tolerance for the cost function, and the cut-off threshold in the hardthresholding process.Tuning the hyperparameters is also vital for the learning outcome, especially the initial value of learning rate α and sparsity constraint λ.The process of hyperparameter tuning is currently done manually with trial and error.
One major limitation of xL-SINDy is the difficulty in designing the library.Prior knowledge of the systems is essential to deciding what candidate functions we should include in the library.A large number of candidate functions in the library are more likely to be sufficient, but it makes the sparse optimization more challenging and less robust against noise [12].Hence, balancing this trade-off is crucial for the outcome of the learning process.
A better mechanism to handle a large library is crucial in applying xL-SINDy to more complex systems with a higher degree of freedom.One possible way to tackle a large number of the library is by the library-bootstrapping method as in the case of Ensemble-SINDy [19].Many smaller libraries are created by sampling the terms without replacement from the original library, and several different models are learned separately.Once the learning process is done, all terms with low probability inclusion are then removed.This process can be repeated until sufficient results are obtained.
So far, we have not considered the presence of external non-conservative force acting on the systems.In real-world scenarios, no matter how small, non-conservative forces such as damping or friction are always present.Taking into account this external force in the model is of importance to make xL-SINDy applicable to real systems.One possible way to incorporate non-conservative force is by using the generalized Rayleigh's dissipation function [20].Like the Lagrangian, Rayleigh's dissipation function is a single scalar quantity and can be incorporated into Euler-Lagrange's equation.We can model the generalized Rayleigh's dissipation function as a linear combination of candidate functions and learn both Lagrangian and Rayleigh's dissipation function simultaneously.

V. CONCLUSION
In this work, we developed a method called extended Lagrangian-SINDy (xL-SINDy) that can discover the true form of Lagrangian of nonlinear dynamical systems from noisy measurement data.We model the Lagrangian as a linear combination of nonlinear candidate functions and use Euler-Lagrange's equation to formulate the objective cost function.We use the proximal gradient method to optimize the cost function and obtain sparse expression of Lagrangian.We demonstrated the effectiveness of xL-SINDy and showed that xL-SINDy is more robust against noise compared to other methods.
It is worth noting that our proposed method out-performs SINDy-PI (parallel, implicit), a recent robust variant of SINDy developed for implicit dynamics and rational nonlinearities.The robustness against noise was improved by 8-20 times compared to SINDy-PI.We believe that xL-SINDy is a promising approach for the identification of interpretable models of nonlinear dynamics.The focus of our next work is to consider non-conservative forces in the model and apply xL-SINDy to real systems.

Fig. 1 :
Fig.1: Block diagram of the proposed method (xL-SINDy).Depending on the case of the problem, a different cost function is constructed.Once the cost function is defined, the cost function is minimized by using the proximal gradient descent method.

Fig. 3 :
Fig. 3: Comparison against true model when external excitation is provided.Training data consists of 100 initial conditions in a time period of 5 seconds each.Validation (extrapolation beyond the training data set) is conducted for 5 seconds afterward.The results shown are taken randomly from one of the initial conditions from the training data set for cart pendulum, double pendulum, and spherical pendulum.

Fig. 4 :
Fig. 4: Simulation results against three different noise levels for three different models: the true model, model discovered by the proposed method, and model discovered by SINDy-PI.Training data consists of 100 initial conditions in a time period of 5 seconds each.Validation (extrapolation beyond the training data set) is conducted for 5 seconds afterward.The results shown are taken randomly from one of the initial conditions from the training data set for cart pendulum, double pendulum, and spherical pendulum.

Fig. 5 :
Fig. 5: Comparison of the maximum level of manageable noise in training data before an incorrect model structure is discovered for three methods: xL-SINDy, SINDy-PI, Lagrangian-SINDy.xL-SINDy provides the most robust performance against noisy training data in all simulations of three different dynamical systems.

TABLE I :
Extracted Lagrangian from simulation data with various noise levels when external inputs are provided.All results are obtained with the computation described in case I.
*Terms highlighted with red color are extra terms that are not supposed to be included in the Lagrangian.All numbers are rounded to 3 decimal places.
Terms highlighted with red color are extra terms that are not supposed to be included in the Lagrangian.All numbers are rounded to 3 decimal places.

TABLE III :
Number of terms used in the library to obtain the model in different dynamical systems.