{'fraction_non_alphanumeric': 0.05170954637779282, 'fraction_numerical': 0.014091062965470549, 'mean_word_length': 4.422900413033502, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Physically Consistent Learning of Conservative Lagrangian Systems with Gaussian Processes

This paper proposes a physically consistent Gaussian Process (GP) enabling the identification of uncertain Lagrangian systems. The function space is tailored according to the energy components of the Lagrangian and the differential equation structure, analytically guaranteeing physical and mathematical properties such as energy conservation and quadratic form. The novel formulation of Cholesky decomposed matrix kernels allow the probabilistic preservation of positive definiteness. Only differential input-to-output measurements of the function map are required while Gaussian noise is permitted in torques, velocities, and accelerations. We demonstrate the effectiveness of the approach in numerical simulation.


I. INTRODUCTION
Euler-Lagrange (EL) systems represent a wide variety of physical systems and are the foundation of a large class of theoretical concepts in control, such as stability, passivity, and in feedback law design. The application of Gaussian Processes (GPs) in learning-based control holds great promise for the control of uncertain systems, improving important aspects like performance [1], safety [2], and dataefficiency [3]. Yet, in general GP regression does not account for physically consistency, hampering the provision of safety guarantees when used in model-based control approaches.
We address this issue by formulating GPs consistent with the relevant physical laws and mathematical models, thus implicitly performing a trade-off between the flexibility of learning and physics-imposed symmetries. Traditionally, parametric techniques [4] have been applied to the identification of EL systems, as reviewed in [5]. However, mostly assuming linearity in the parameters, these methods eventually become infeasible with increasingly nonlinear, and coupled, dynamics, as for example in soft robots [6], aircrafts [7], and marine vehicles [8].
Based on Bayesian principles, GPs allow the inclusion of analytical prior knowledge, thus forming a bridge to parametric model-based approaches. The authors of [9] propose the regression of Hamiltonian systems with GPs, preserving the symplectic structure, but without any guarantees and requiring the availability of noiseless impulse measurements. A Lagrangian kernel is derived in [10] within the framework of Reproducing Kernel Hilbert Spaces (RKHS). While analogies to our stochastic GP-based approach naturally exist, the work lacks a further specification of the associated function space, prohibiting a comprehensive analysis of control relevant properties closely intertwined with physical consistency. In [11], the EL residual dynamics are compensated by employing independent GPs per degree of freedom, neglecting correlations and structure. This leads to an implicit differential equation system, the necessity of noiseless acceleration measurements and the curse of dimensionality. Guaranteeing physical consistency of data-driven methods still represents a major open problem.
The main contribution of this work is the derivation of a physically consistent Lagrangian-Gaussian Process (L-GP) constructed by linear operators, embedding the differential equation structure into the model. Combined with quadratic kernels, we extend the GP concept from scalar function to symmetric matrix spaces. Based on a Cholesky decomposition of the underlying matrix kernel, we are able to provide a probabilistic guarantee for the positive definiteness of the estimates. Further physical consistencies of the L-GP analytically preserving equilibria, passivity, and conservatism are formally shown.
After defining the problem setting in Section II, we review the framework of GPs in Section III and describe our proposed unification with EL systems in Section IV. Section V presents our main results. Numerical illustrations validate theory and effectiveness of the approach in Section VI.

II. PROBLEM FORMULATION
In this paper, we consider the class of conservative EL systems whose equations of motion are given by 1 based on the set of generalized coordinates q ∈ R N and the unknown Lagrangian L. The generalized forces τ ∈ R N are the measurable, possibly underactuated, inputs to the system. Assumption 1: The unknown Lagrangian is autonomous with L ≡ L(q,q) = T (q,q) − V (q) and composed of the 1 Notation: Vectors a and matrices A are denoted with bold lower and upper case characters, respectively. The nabla operator ∇ multiplied with a differentiable scalar field f (x) gives its gradient w.r.t. x ∈ X . The partial gradient w.r.t. a subset of variables y ∈ Y, Y ⊂ X , is described by ∇yf . I denotes the identity, 0 the zero and 1 the ones matrix. A collection of D indexed vectors x i ∈ R N is stacked row-wise in the matrix ∈ R D×N such that x ij refers to the j-th component of x i . For notational convenience, we omit the indexing in the matrix expressions and Cov[·] denote the expectation, variance and covariance operators w.r.t. random variables. The Hadamard product and square are denoted by x • y and x •2 = x • x, respectively, the Kronecker product with ⊗, and the direct sum with ⊕. With λ k (·) we denote the k-th largest eigenvalue corresponding to the components of the decreasingly ordered eigenvalue vector λ ↓ (·). difference between unknown kinetic energy T : R N × R N → R and unknown potential energy V : R N → R.
Assumption 2: The kinetic energy is of the form T (q,q) = 1 2q T M (q)q where M : R N → R N ×N is the unknown, (symmetric) positive definite, inertia matrix.
Assumption 3: The potential energy is decomposable into the sum V (q) = G(q)+U (q) of gravitational and elastic potential energies G and U , respectively, where the former has an equilibrium, G(0) = 0, ∇ q G(0) = 0, and the latter is given by U (q) = 1 2 q T S(q)q with unknown, (symmetric) positive definite, stiffness matrix S : Assumptions 1-3 essentially describe the considered system class and represent its physical properties. Assumption 1 restricts the systems described by (1) to be only implicitly dependent on time and free from disturbances caused by other external forces such as dissipation. Extending to handle dissipative and time-variant systems is possible and topic of current research. Assumption 2 is a valid approach when dealing, e.g., with classical mechanical rigid body systems, while Assumption 3 constrains chosen coordinates to have an equilibrium at the origin.
Exploiting Assumptions 1-3 and applying the chain rule to (1), we obtain the well-known matrix-vector expression introducing the (non-unique) generalized Coriolis matrix C(q,q) : [12] and the vector of generalized potential forces g(q) Note that there are other ways to define the matrix C(q,q). This particular choice, however, has the unique property thatṀ − 2C is skew-symmetric.
Our objective is to approximate the unknown Lagrangian function L(q,q) by a data-based estimateL(q,q) which is physically consistent, i.e., maximally aligned with the existing physical knowledge, in our case Assumptions 1-3. We assume access to noisefree positional observations q i , i = 1, . . . , D with D ∈ N, and to potentially noisy velocity, acceleration and torque measurementṡ Remark 1: Extending to noisy positional inputs is possible via certain GP-based methods as in, e.g., [13], [14], yet they are mostly based on first-order Taylor approximations.

III. GAUSSIAN PROCESS FRAMEWORK
This section gives a brief overview of the mathematical framework of GPs based on [15]. For a more complete introduction, the reader is referred to the literature [16]- [18].

A. Inference with GPs
For x, x ′ ∈ X in a continuous domain X ⊆ R 2N , a GP with mean m(x) and covariance k(x, x ′ ), denoted by and with scalar f (x) ∈ R, is a stochastic process extending the Gaussian probability distribution from random variables to functions. As such, it inherits the convenient mathematical properties of the Normal distribution. In particular, any finite marginalization z = [z i ] of function observations z i = f (x i ) + ζ i , corrupted by the zero-mean white-noise process ζ i ∈ R and corresponding to the finite subset of function evaluations at {x i } ⊂ X , is jointly Gaussian distributed. A GP is fully characterized by its mean and kernel function m( respectively. The former allows the analytical inclusion of prior knowledge about the unknown function, whereas the latter determines higher level functional properties such as smoothness. Estimation of the covariance describing hyperparameters is mostly done via optimization of the marginal likelihood to maximize the probability of observing the measured outputs. For regression, GPs exploit the joint Gaussian distribution of observations z at measurement locations X = [x ⊤ i ] and prediction f (x) w.r.t. the test input x ∈ X given by where we have introduced the mean m = [m(x i )] and covariance k( Here, we have dropped the dependencies of m(X), k(X, x) and K(X, X) on X for notational simplicity. Conditioning on the observations yields a conditional Gaussian distribution with posterior given analytically bŷ

B. Linear Operators and GPs
Derivatives and integrals of GPs remain GPs, since these operations are linear. Applying a linear transformation operator T x to the GP of (4) leads to a new, possibly multidimensional, GP in the form of with transformed mean m T (x) = T x m(x) and kernel Here, T ⊤ x ′ is applied from the right [9] as is necessary when considering matrix operators.

IV. REGRESSION OF LAGRANGIAN SYSTEMS
Our key idea is to consistently model the unknown Lagrangian using a GP. We propose a structural approach for the application of GP regression (4)-(6) to EL systems (1)-(2): deriving specific GPs and kernels, we ensure consistency with physical requirements, as in Assumptions 1-3, and embed the structure of the differential equations using Nabla operators ∇, thereby inducing multidimensional L-GPs according to (7)-(8) based on the linearity of differentiation.

A. Energy Structuring
Firstly, we model the Lagrangian as a composite GP with the subjacent kinetic and potential energy GPs Here, the means m L = m T −m V and covariances k L = k T + k V specify the L-GP (9) in a structurally consistent manner.
In particular, the potential energy GP (10b) depends only on the position q by construction. It is split up further into the sum V (q) = G(q) + U (q) of independent gravitational and elastic potential GPs (11b)

B. Quadratic and Positive Definite Structures
Consistency of kinetic and elastic GPs is ensured by according design of their mean and covariance functions. If available, we impose on the according priors the same quadratic structure based on, e.g., for the kinetic energy, the (symmetric) positive definite a-priori matrix M 0 (q) ≻ 0, ∀q ∈ R N , as described in the following.
Assumption 4: The priors m f k ∈ {m U , m T } are in the same form as the unknown energy functions In addition, the gravitational prior fulfills the equilibrium condition such that m G (0) = 0 and ∇ q m G (0) = 0 hold. Having specified the means in (10a) and (11b), we now move on to the kernels. Consider the quadratic functional with differential-operational index k ∈ {0, 1} and the (symmetric) positive definite Cholesky decomposed matrix kernel The Cholesky factor R k is an upper-right triangular matrix consisting for 0 ≤ n ≤ m ≤ N of kernels r knm (q, q ′ ). We now use (13) to set the covariances in (10a) and (11b) to: Thus, we implicitly model the stiffness and inertia matrices denoted in accordance with (12) by H k := [h knm ] as symmetric, matrix-valued, GPs with entries where h knm = h kmn , and analagously constructing H k0 := [η knm ]. The kernels θ knm = ρ ⊤ kn ρ km stem from the inner product of the non-zero columns from (15), i.e, ρ km := [r klm ] ∈ R P for P = min{n, m}. Note that, due to the nonlinearity of multiplication, the r klm may not be associated any further with underlying GPs, or more specifically, Gaussian Random Variables (GRVs). However, the covariances θ knm remain valid kernel functions according to [15].

C. Differential Equation Embedding
Let us now use the chain rule to expand L q ≡ d dt ∇q −∇ q from (1) by writing Applying the Lagrangian-differential vector operator (17) to the GP from (9), we obtain for the torques (1) with mean vector m τ = L q m L and kernel matrix K τ = L q L ⊤ q ′ k L . As indicated in (18) via the semicolon,q andq have the role of regressors, whereas q remains a conventional input.
dropping dependencies on all variables for notational brevity. Also, we have introduced: the stacked vector of outputs y = vec(Y ⊤ ) with mean m y = [m τ (q i ;q i ,q i )] and covariance the energy variances σ 2 T (q,q) = k T (q,q, q,q) and σ 2 V (q) = k V (q, q), the equilibrium variance and covariance K y0 = [k yV (0), ∇ ⊤ q k yV (0)], as well as the remaining equilibrium-potential k 0V (q) and Lagrangiandifferential energy k yT (q,q), k yV (q) covariances The fully formulated L-GP (19) is the first part of our contribution. Before moving on to analyze its physical consistency, we propose a method for partial input noise compensation.

E. Nonlinear Noise Compensation
The regressor structure stemming from the application of (17) to (16b) allows for the compensation of noise in the differential inputsq i andq i . Therefore, as a last step, we combine (2) and (18) in order to transform (3b) into Approximating the product of two GRVs to remain Gaussian, we obtain a heteroscedastic (state-dependent) noise model where with the inertial variance Θ 1ii := Θ 1 (q i , q i ), the Coriolis i ](∇ •2 qi ) ⊤ matrices as well as the main diagonals vectors σ αi and σ ωi of Σ αi and Σ ωi , respectively.
Proposition 1: The distribution (22) with state-dependent variances (23) is a stochastically consistent approximation which becomes exact for the limit case Σ αi , Σ ωi → 0, or the case of certainly known T = m T with k T ≡ 0.

V. INVARIANCE PROPERTIES
Having introduced the novel framework of L-GPs, we now investigate their properties and derive certain guarantees.

A. Quadratic Form
Employing the kernel structure (13), quadratic form of the energy GPs (10a) and (11b) can be guaranteed analytically. In order to specify our result in Lemma 1, we require the following definition: for δx i ∈ R N and fixed q, the directional derivative of Θ k (q i , q) along δx i at q i gives The matrix Φ k is the directional matrix-derivative of Θ k from (14) moving through position q i with velocity δx i . In particular, note that holds for fixed q by deriving w.r.t. the virtual time τ . Lemma 1: The mean estimatesf k (q) ≡ E[f k (q)|y, γ 0 ] of kinetic and elastic GPs (10a) and (11b) in the joint model (19) are analytically guaranteed to be of quadratic form: The posterior matrix estimatesĤ k are decomposed aŝ where the basis matrices N ki are given by with the subvectors δx i = [∆x n ] iN n=1+(i−1)N of the transformed innovation ∆x = K −1 D ∆y, the Kalman-like kernel gain matrix K D = K y − K y0 K −1 0 K ⊤ y0 , and the innovation difference ∆y = y − m y .
Proof: We set G ≡ 0 w.l.o.g., leading to V = U ≡ f 0 and k 0V = 0 due to (16a) and (13). Marginalizing (19) Then, applying standard inversion formulas for partitioned matrices together with γ 0 = 0, we can writê as the application of (17) to (13) such that k yf k = [l ki ], it can be verified that holds after some differential vector-algebraic computations, where we utilize among other intermediate steps that The full derivation is omitted due to space limitations. Further exploiting the commutativity of the Hadamard product and the definition of the projected gradient matrix Φ k in (24), we follow that δx ⊤ i l ki (q) = 1 2 q (k)⊤ N ki (q)q (k) with N ki as in (27). Then, eliminating the skew-symmetric part of N ki combined with the quadratic priors m f k from Section IV-B, we reformulate (28) using k ⊤ yf k K −1 D ∆y = i δx ⊤ i l ki and obtain (25)-(26).
Essentially, Lemma 1 follows from the quadratic kernel structure (13). The summation (26) represents an extension of matrix decompositions from vector to matrix spaces which are aligned with the Lagrangian structure (1)-(2).

B. Positive Definiteness
We now move on to investigate positive definiteness of the matrix estimatesĤ k . Before that, however, we impose a minor structural restriction on the used kernel functions.
Assumption 5: The Cholesky matrix kernel components r knm from (15) and the gravitational covariance k G from (11a) belong to the class of metric kernels M given by parametrized by σ knm , Λ knm and σ G , Λ G , respectively. Note that class M in Assumption 5 includes a wide variety of covariance functions including the commonly used squared exponential kernel [15]. Therein, kernels k(q, q ′ ) ≡ k M (|d|) are Gaussian radial basis functions. They depend only on the squared length of the difference in the Riemann space [20, p. 243] defined by metric Λ ≻ 0. The special case of variance We can now explicitly state the components r knm of the Cholesky factor R k from (15). Introducing for 0 ≤ n ≤ m ≤ N the functionsr knm ∈ M from Assumption 5 with identical Riemannian hypermetrics Λ knm = Λ k , we set We now make a stochastic analysis of the data set D gathered via (3a)-(3b), enabling the probabilistic investigation of positive definiteness preservation. Based on the product p we define for µ k , ν k ∈ R D the composite random variable for those i ∈ {1, . . . , D} which satisfy component-wise ∀n ∈ {k, 2k}, k ∈ {0, 1}, the condition p (n) i > 0, by Here, the projection ϕ ki , the differential angleθ ki and the normalizing function Υ k (p i ) are given by whereĤ ki = 1 2 (N ki +N ⊤ ki ) from (26). The random variable β k from (33)-(35) is the result of nonlinearly transformed GRVs stemming from the data set D. We now formulate the main result.
Theorem 1: Consider the posterior, generalized, matrixvalued L-GP estimateĤ k (q) from (26) with metric covariances r knm ∈ M forming the Cholesky kernel Θ k (q i , q) from (14) identically specified by Riemannian hypermetric Σ −1 d k via Assumption 5 and (31). Positive definiteness of the generalized matrix estimate is guaranteed with probability where the physical inconsistency measure δ k (q) is upper bounded according to which depends only on the absolute values of the positional joint distances |d j (q)| = |q − q j |.
Proof: Firstly, we express the basis matrices (27) from Lemma 1 writing Θ k ≡Θ k (d) based on Assumption 5 and (31). Proceeding analogously via Φ k ≡Φ k (d, δx i ), we follow thatΦ with the projection ϕ ki and the difference vector d i as in (34a) and (32b), respectively. Then, we can write for each matrix summand N ki in (27) that dropping dependencies on q for the sake of notational brevity. Consider now the i-th summandĤ ki = 1 2 (N ki + N ⊤ ki ) from (26). Exploiting where Σ 2 f k = Σ ⊤ f k Σ f k , and its eigenvalue decomposition we now leverage the eigenvalue behavior of the sum of Hermitian matrices [21] in combination with Weyl's inequality [22,Theorem III.2.1]. Thus, with the Hermitian summation H k = iĤ ki as in (26), we derive Then, utilizing the product constructions (37), we obtain by successively exploiting the invariance of eigenvalues w.r.t. similarity transforms and by using d i as well as the definition of the random vector β k ∈ R D as given component-wise by (33)-(35). Taking the complement of the set in which the lower bound (38) is positive, we arrive at (36) and finish the proof.
Note that, for a given data set D, the upper bound on δ k in Theorem 1 can be used to deterministically estimate the domain in which positive definiteness is fulfilled. A probabilistic computation can be made by means of the distribution of products and quotients of continuous, nonstandardized, independent Normal random variables [23].

C. Equilibria
Despite the point-wise inclusion of γ 0 = 0 in the distribution (19) being stochastic, we can provide an analytical guarantee for the preservation of the equilibrium at the origin.

D. Energy Conservation
Lastly, the L-GP is also equivalent to a dynamic system. We investigate this intuition in the following.
We now show that the equivalent L-GP dynamics (42) preserve the physical property of energy conservation as a last result. For notational consistency, we use the (estimative) state vectorx ⊤ = [q ⊤ ,q ⊤ ] in the following.

VI. NUMERICAL ILLUSTRATIONS
In this section, we now demonstrate the efficacy of our approach and validate the derived theoretical results, choosing a simple example for the sake of easy interpretation.

A. Setup
We benchmark our method on the two-link robotic manipulator from [12, p. 164]. Gravity g = 10 ms −2 acts along the positive x-axis, leading to an equilibrium at the origin q = 0. The links have unit masses m n = 1 kg and unit lengths l n = 1 m for n ∈ {1, 2}. Estimates are available but erroneous:m n = (1 + ∆ n )m n ,l n = (1 + ∆ n )l n , where ∆ n = (−1) n /2.

B. Results
As a first experiment, we demonstrate the closed-loop applicability of the L-GP model. For this, we compare a standard PD tracking controller [12, p. 194] with its L-GP-based version as shown in Fig. 1. Starting from the equilibrium, the controllers with gains K p = K d = 10I have the task of following a sinusoidal reference trajectory q d (t) = a d sin t1 with amplitude a d = π/2. The simulation shows that the L-GP-based version considerably improves the accuracy, demonstrating reliable performance despite the suboptimal distribution of the data points, independent from the reference trajectory, requiring substantial extrapolation.  Next, we validate the physical consistency of the L-GP. Therefore, using the same parametrization as in the previous experiment, we evaluate the mass inertia estimate over the joint domain q ∈ [−a d , a d ] 2 and investigate its eigenmanifolds, as shown in Fig. 2. Clearly, the L-GP accurately approximates the positive definite function space despite only being a subcomponent of the input-output relation. The consistency of the gravitational potential estimate with equilibrium is validated in Fig. 3, along with the conservatism of the equivalent L-GP-based dynamics simulated for different initial conditionsx(0) = [a 0 1 ⊤ , 0 ⊤ ] ⊤ with displacement amplitudes a 0 = 0.1 (red line), a 0 = 0.5 (yellow line)

VII. CONCLUSION
This paper presented an approach for the identification of uncertain Lagrangian systems exploiting prior physical knowledge by kernel construction. Physical consistency of the data-driven method in terms of guarantuees for the fulfillment of certain properties was proven rigorously and validated in numerical simulation, along with the effective applicability of the model to an exemplary tracking control problem. Future research will focus on the extension of the approach to handle dissipative and time-variant systems. Another promising direction is the application of the Lagrangian variance estimate for uncertainty quantification.

VIII. ACKNOWLEDGMENTS
This work was supported by the Consolidator Grant "Safe data-driven control for human-centric systems" (CO-MAN) of the European Research Council (ERC) under grant agreement ID 864686, and by the Horizon 2020 Research and Innovation Action project "Rehabilitation based on Hybrid neuroprosthesis" (ReHyb) of the European Union (EU) under grant agreement number 871767. The authors gratefully acknowledge the thoughtful comments from L. Evangelisti, Institute of System Dynamics and Control, German Aerospace Center (DLR), and the fruitful discussions with A. Lederer, Chair of Information-oriented Control, Technical University of Munich.