Sparse Identification of Lagrangian for Nonlinear Dynamical Systems Involving Dissipation Function

More than three hundred years ago, Isaac Newton tried to figure out what made an apple fall off the tree; behind the phenomenon, he found gravity. That is one of our first attempts to find the system dynamics through the phenomenon. Since then, modeling has become one of the eternal topics in science research; even now, scientists are still trying to model the complex relationships between stars. However, there are too many dynamic systems in daily life, and it is impossible to model all of them manually. Still, knowing the dynamics might be important since it helps us to understand how the system's parts contribute to the movements. That makes the autonomous modeling method helpful for saving efforts. However, the noisy, complex environment makes real-world autonomous modeling more challenging than simulation, which is still an open topic [1–4].

There is a debate between complexity and parsimony in the quest for precision in modeling. Complex models offer detailed descriptions through complex equations but are often hard to interpret as explainable knowledge. A representation is the neural network model, which can fit the data with high accuracy, but it is hard to find analytical meaning from the network. Conversely, parsimonious models aim for simplicity and capture the system's essence with minimal components. As a representation, the equation of motion describes the system's movement and provides insights into the system's structure. Despite the difficulty of deriving such models, contemporary research prefers parsimonious models for their analytical elegance [5–11].

Though studies in this area started more than a hundred years ago, applicable methods only began to appear after machine learning tools were developed [12–15]. In 2007, Lipson and co-worker explored merging machine learning into the autonomous modeling process [6]. In 2016, Brunton et al. proposed the game-changing sparse identification of nonlinear dynamics (SINDy) algorithm [7]. The machine learning method that removes the irrelevant terms and gives a sparse result, commonly called sparse regression, is the key that helps SINDy to lower the computation cost. That made SINDy to be the first applicable modeling algorithm. It was evaluated in many areas, from chemical kinetics to chaos theory [7,16–24].

With its proven performance, we choose it as the basement of this research. However, vanilla SINDy is sensitive to noise, so applying SINDy in the real world is challenging. To mitigate this issue, Brunton and co-workers proposed a more robust variant called SINDy with parallel and implicit expression (SINDy-PI) [25]. However, with the maximum capable Gaussian noise level as $10−3$⁠, it was still not robust enough for real-world tasks. To make the algorithm even more robust, our group proposed the extended Lagrangian-SINDy (xL-SINDy) method, which changes the learning target from the equations of motion to Lagrangian, which describes the system dynamics in a simplified way [26,27]. Thanks to the simplicity of Lagrangian, the xL-SINDy can manage noise magnitude up to $10−1$⁠, outperforming SINDy-PI. However, as a tradeoff, xL-SINDy cannot manage those systems with energy exchange with the outer environment, commonly called nonconservative systems. Since a typical energy exchange is the dissipative force and it is everywhere in the real world, most dynamical systems in the real world are nonconservative. We need to extend the xL-SINDy to manage the nonconservative case before we apply it in the real world to solve the autonomous modeling challenge.

To better understand the relationship between SINDy, SINDy-PI, xL-SINDy, and our method, we can disassemble the whole recognition process into three parts:

1. The base algorithm: This determines the computational efficiency of the process.

2. The identification object: The object represents the system dynamics and works as the identification object; it affects the process robustness.

3. Formation of the regression loss: This affects the optimization criteria and brings the nonconservative capability of the process.

All four methods use the same base algorithm—the SINDy algorithm. For SINDy-PI, the identification object is replaced, providing an enhanced ability to handle greater robustness. xL-SINDy improved with the identification object and loss formation, with robustness surpassing SINDy-PI; however, it cannot handle nonconservative cases. In our research, we utilize the xL-SINDy framework but modify the loss to improve its capability, thereby retaining xL-SINDy's superior robustness while overcoming the capability limitations. Figure 1 illustrates the comparison between SINDy and its variants.

In this paper, we combine a nonconservative expression called the Rayleigh dissipation function and the xL-SINDy, providing a highly robust and nonconservative capable solution for system identification and one step closer to solving real-world problems.

This section introduces the vanilla SINDy and xL-SINDy and then discusses the modification to adopt a nonconservative system. A training efficiency test helps select the appropriate algorithm for the sparse regression task. Finally, the workflow of the whole algorithm is introduced, and a figure demonstration is given.

As the basement of all other SINDy variants, vanilla SINDy enables efficient system dynamics modeling with a sparse regression procedure. Figure 2 explains the algorithm step-by-step. The simplest form separates the equations into the coefficients and system states' functions. Then, the functions are stored in a vector, and the coefficients are stored in a matrix; the multiplication of these two can regenerate the equations. The next step assumes that the coefficients are unknown but that data from the system are given. In this case, the coefficient matrix would be initialized with random values, and a linear regression could help find the coefficients from the data. The regression loss is the error between the acceleration prediction and the true value. Finally, for the complete SINDy case, the coefficients and system states' functions are unknown. The initial function vector, also called the candidate library in SINDy, collects all possible combinations of the system states. The coefficient matrix would get randomly initialized, similar to the second step. The algorithm uses sparse regression to discover the relevant terms inside the huge function vector. It promotes most of the coefficients to zero, finding the dynamics from the data.

The regression process is to find a nonzero solution $C′$ to make this equation stand, and its multiplication with the $Φ$ gives the final learned Lagrangian. The defect of xL-SINDy comes from the zero on the right-hand side of Eq. (5), limiting it to only describe conservative systems.

To address the limitations of the previously described xL-SINDy method, which pertains specifically to conservative systems, we extend the capability of the Euler–Lagrange equation by adding a nonconservative term on the right-hand side. Specifically, we introduce the Rayleigh dissipation function model, which provides a universal solution for modeling various common types of dissipative force. It can generate different types of dissipative force by changing only one hyperparameter, being able to produce a large number of candidate functions in advance; it aligns with SINDy's need for the candidate function vector construction. By adding the Rayleigh dissipation function terms into the candidate function vector, our algorithm can learn the system's Lagrangian and its nonconservative parts at the same time. Although the Rayleigh dissipation function already contains the common types of dissipative force, it is necessary to admit that the Rayleigh dissipation function cannot cover all situations. Adding additional terms to the candidate function vector can address this issue.

The left side of this equation is the same as Eq. (3), with the Rayleigh dissipation function term added on the right-hand side. With a minus term as the system energy flow, the system loses energy due to dissipative force, which is the nonconservative case we want to describe.

Assuming we are working with a double pendulum system, $q˙j$ refers to the velocity of the $jth$ joint, and $cj$ refers to the coefficient of the $jth$ joint's dissipation function. The n refers to the order of the dissipation function. With such a function, we can implement different kinds of dissipative force into the system by changing n. For example, when $n=0$⁠, the function refers to a constant dissipation function like constant dissipative force, and when $n=1$⁠, it refers to a dissipation function that is linear to one of the system states like linear dissipative force. All different components can be easily added up to compose the whole dissipation function. In the experiment, the candidate function vector includes dissipative force created by the Rayleigh dissipation function for every system state with n ranging from 0 to 3, while a constant or linear type dissipative force is applied in the system.

After understanding how the Rayleigh dissipation function works, the next step is to add it to the xL-SINDy framework. For convenience in sparse regression, we need to divide the Lagrangian and the Rayleigh dissipation function into the candidates part and coefficient part as follows: $L=∑k=1pckϕk,D=∑j=1hejψj$⁠.

Here, the $ϕk$ and $ψj$ refer to the Lagrangian and dissipation function candidates, respectively, and $ck$ and $ej$ refer to the coefficients of these two candidates' libraries.

This demonstration mainly focuses on the passive scenario in which no active force is applied to the system, and no prior knowledge is assumed. With some modification, this method can also be applied to active systems and systems with partial prior knowledge. An active model can be obtained by including the dissipation function as part of the external force. The model with partial prior knowledge can be obtained by separating the dissipation function on the left-hand side of the equation.

The $τext$ refers to the known external force input to the system.

Having some prior knowledge means knowing certain candidates in the expression and their coefficients. Here, $Ml,Nl,Ol$ refers to the components formed by the known candidates.

The sparse regression training in the SINDy algorithm is challenging for its high sparsity compared to other common applications of sparse regression algorithms. For standard sparse regression cases, the ratio of irrelevant functions inside the whole vector, commonly named sparsity, is around 80% to 50%. In SINDy tasks, the sparsity is usually higher than 90%; such a sparse dataset is challenging to fit. To address this challenge, a robust, efficient, sparse regression algorithm is needed. In experiments, the proximal gradient descent (PGD) [28] methods perform well in high-sparsity sparse regression thanks to their robustness brought by the proximal operator. To find the optimal option, we tested three variants of PGD and vanilla PGD; the result is shown in Fig. 3. Although the PGD with line-search (PGD with Ls) [29] seems to converge with the least iterations, it takes longer for each iteration than all other methods. As a result, its actual time consumption sits between the accelerated proximal gradient descent (APGD) [30] and fast iterative shrinkage threshold algorithm (FISTA) [31], making the APGD the best choice for our task.

In this section, we tested the algorithm's robustness on four different systems and then verified its generalizability with two extra tests. With a maximum capable noise level as $10−1$⁠, the proposed method maintains the same robustness level as the xL-SINDy. In the generalizability test, the algorithm shows an accuracy of over 90% across different situations, proving its potential as a universal solution.

For ease of reading, the paper only compares the simulation curve between the true and obtained models. A table of all obtained expressions is in Table 1. In the single pendulum, $θ$ refers to the angle of the pendulum; in cart pendulum, x refers to the displacement of the cart and $θ$ refers to the angle of the pendulum; in the double pendulum, $θ1$ refers to the angle of the first joint which is the joint that got fixed with the environment, $θ2$ refers to the angle of the second joint which is the joint that connects the two bars; in the spherical pendulum, $θ$ refers to the polar angle, and $ϕ$ refers to the azimuthal angle.

We modified the training process with two special mechanisms besides the mentioned measures: two-phase training and scaling [32]. This accelerates and stabilizes the training process and does not affect robustness. With a weak connection with the main algorithm and page limitation, we skip the details of these mechanisms.

With Fig. 4 summarizing the whole algorithm, here is the training flow: First, select the Lagrangian and Rayleigh candidate functions, then randomly initialize the corresponding coefficients. By applying them to the Euler–Lagrange equation, we obtain the system's energy flow expression. Next, substitute the system states into Eq. (13), (14), or (15) to predict the system's energy flow. Using the error between the prediction and the true value, update the coefficients and fit them to the data, which gives the final Lagrangian and Rayleigh dissipation functions. We also provide some visualizations of the libraries in Figs. 10–13; these figures are provided in the  Appendix.

In the results, we noticed that the coefficients for the dissipative function tend to decline with the growing noise magnitude, a possible explanation would be the dissipative force' effect is gradually being covered with the growing noise, making it harder to recognize the dissipative part under high noise level.

To evaluate the effectiveness of our method, we selected single pendulum, cart pendulum, double pendulum, and spherical pendulum systems as subjects for our experimentation; the detailed configurations of these systems are defined in our previous work [27] and shown in Fig. 5. These choices were made due to their relatively simple but chaotic behavior, making them ideal candidates for comparison with SINDy-PI. The dynamics governing these systems are defined by their respective Lagrangian, as presented in Table 1. To include dissipative part in the system, a linear dissipative force is applied to all four systems, as the dissipation function shows in Table 1. The whole process is based on simulation, and the data sample frequency is 50 Hz.

To evaluate the robustness of the proposed method, we added Gaussian noise into the generated data. We aligned our noise magnitude settings with those of xL-SINDy for comparison convenience, selecting values of $10−4$⁠, $10−3$⁠, $2×10−2$⁠, $6×10−2$⁠, and $10−1$⁠. The tests reveal that the proposed method maintains the same robustness level as the xL-SINDy and outperforms SINDy-PI but with the extra capability to model nonconservative systems.

To show that the algorithm can perform stably with different pendulum systems, we also tested the algorithm under the same model type but with different physical parameters and dissipative force types.

For the sake of readability and to aid in interpreting results, all Lagrangians in this study have been scaled relative to the first coefficient. This approach is valid, as scaling the Lagrangian by an arbitrary constant does not alter the system dynamics. For the same reason, systems under different noise levels are started at the same initial position in the figures below, though in real experiments the position is randomly initialized. In the figures, the first 5 s would be the training data with noise and the obtained model's data, while the rest would be the validation data that started from a different initial position and the obtained model's data. By showing both the training and validation comparison to the obtained model, it is proved that it is indeed a universal model that fully describes the system's dynamics.

For the single pendulum system, we attempt to recover the system dynamics with the Rayleigh dissipation function model. In sample generation, the pendulum is randomly initialized in the range $[−π,π]$ for the angle and zero for the velocity. The candidate function library includes polynomials and trigonometric functions up to second order for the Lagrangian, and dissipation functions for n ranging from 0 to 3, eventually building a library containing 17 functions, with three true functions inside the library (two from the Lagrangian and one from the dissipation function); the task sparsity is 17.64%. Figure 6 gives an intuitive simulation based on the result of the Rayleigh dissipation function incorporated method. The robustness of the algorithm is also proven. Although not shown here, for the conservative system, the coefficient shifting is less prominent along with the growing noise; the same magnitude of noise has a greater influence on nonconservative system. This effect stems from a larger relative error when the oscillation magnitude decreases at the late stage of nonconservative systems. Overall, the positive result in this simple model proves the algorithm's viability, giving confidence for tests on complex systems.

The angle $θ$ of the pendulum is randomly initialized at the range $[−π,π]$⁠, and the cart position is zero, both with zero velocity. The candidate function library includes polynomials and trigonometric functions up to fourth order for the Lagrangian, and dissipation functions with n ranging from 0 to 3, eventually building a library containing 61 functions, with six true functions inside the library (four from the Lagrangian and two from the dissipation function); the task sparsity is 9.84%. The simulation result is shown in Fig. 7. As the results show, our method can find the correct candidates for noise levels up to $10−1$⁠, although with some coefficients shifting, the simulation curve of our obtained can still follow the true model.

We also performed a comparative analysis using the well-known SINDy-PI method to validate our proposed algorithm further. In Fig. 7, the dotted lines give the results of SINDy-PI. As the results indicate, SINDy-PI method is more sensitive to noise. Even though it is claimed that SINDy-PI can handle noise magnitudes up to $10−3$⁠, the actual robustness of the algorithm is also affected by the system complexity. The involvement of dissipative force significantly adds complexity, causing SINDy-PI's performance degradation. In our cart pendulum test, SINDy-PI can handle only up to $10−4$ noise magnitude. In contrast, our proposed method can accommodate noise magnitudes up to $10−1$⁠, matching the robustness of the original xL-SINDy algorithm.

The pendulum is randomly initialized for sampling, both joints' angles are in the range $[−π,π]$⁠, and velocities are zero. The candidate function library has 95 candidates, with eight true functions in the library (six for the Lagrangian and two for the dissipation function); the task sparsity is 8.42%. Testing with the double pendulum system in Table 1, Fig. 8 offers an intuitive simulation-based comparison between the model generated by our algorithm, the model generated by SINDy-PI, and the true model. Our algorithm effectively recovers the correct candidate terms up to a noise magnitude of $10−1$⁠, though there is a noticeable shift in coefficients when the noise level reaches $6×10−2$⁠. This could be attributed to the relative noise level issue discussed in Sec. 3.2. Still, the algorithm exhibits $103$ times higher robustness than SINDy-PI.

Like the cart pendulum case, SINDy-PI could only handle noise level up to $10−4$⁠. For higher noise level, the SINDy-PI can only discover part of the true model or even fail to capture any essence of it. In the meantime, our method can still find the correct candidates and follow the true model's curve at the noise level of $10−1$⁠, further proving the robustness of our method.

Finally, we also tested our method using a spherical pendulum mentioned in Table 1. With a polar coordinate definition, the pendulum's polar angle $θ$ is initialized in the range $[π/3,π]$⁠, while the initial azimuthal angle is zero, and the angle velocities are both zero. The candidate function library has 65, with five true functions in the library (three for the Lagrangian and two for the dissipation function); the task sparsity is 7.69%.

The test results are shown in Fig. 9; compared with the other three cases, the spherical pendulum turns out to be more challenging to recognize. Our method can only handle noise levels up to $6×10−2$⁠, as SINDy-PI could not give any valid result even for cases without noise.

To evaluate the algorithm's generalizability, the algorithm is also tested under varying physical parameters. To be more specific, we utilize the double pendulum system and change its mass, length, and damping coefficients. Table 2 provides a visual representation of the model's accuracy through a heatmap under a noise level of $10−3$⁠; the deeper color refers to a higher accuracy. The accuracy is calculated using the mean square error of the coefficients between the identified model and the true model. Here, $l1,l2$ refer to the length of the pendulums, $m1,m2$ refer to the masses, and $b1,b2$ refer to the dissipative coefficients. Across various physical parameters, our algorithm consistently recovered system dynamics with an accuracy exceeding 90%, affirming its promising applicability in diverse settings. Besides, model performance under different dissipative force types and magnitudes is also tested. The test contains constant and viscous dissipative force, with coefficients ranging from 0.05 to 0.5. With results shown in Table 3, the overall accuracy stays above 90%, showing the algorithm's capability. The result also shows that the algorithm performs better under constant dissipative force. Overall, the proposed method generalizes different system settings and dissipative force, proving its potential as a universal solution.