Structure-preserving formulations for data-driven analysis of coupled multi-physics systems

We develop a novel methodology for data-driven simulation of coupled multi-physics systems. The result of the method is a learned numerical integrator of the coupled system dynamics. In order to preserve the fundamental physics of the coupled systems, and thus preserve the geometrical properties of the governing equations—even if they may be completely unknown—we impose a port-metriplectic structure on the system evolution, i.e., a combination of a symplectic evolution for the system energy with a gradient flow for the entropy of each system, which can be exchanged through predefined ports. The resulting method guarantees by construction the satisfaction of the laws of thermodynamics for open systems, leading to accurate predictions of the future states of their dynamics. Examples are given for systems of varying complexity, based on synthetic as well as experimental data.


Introduction
The analysis of complex physical systems from experimental data is a highly topical subject with countless practical applications.Among them, we cite the development of digital twins [1,2] or structural health monitoring [3], to name but a few.Machine learning techniques developed in recent years shed light on the behavior of these systems, particularly when there is no detailed knowledge of the physical laws governing their behavior.
Although great advances have been made in the development of these techniques, there are still certain difficulties that have prevented their widespread use in industry.These include those related to the accuracy and stability of the predictions made, which are often sensitive to variations in the input data.
This has led to interest in the development of techniques that can ensure a certain degree of accuracy, as well as compliance with known physical laws.If these are fully known, in the form of PDEs, the most widespread technique is the use of physics-informed neural networks, PINNs [4].However, there is often a discrepancy between the predictions made by these physical laws and the actual behavior of the system under analysis [5].On other occasions, these equations are not known and these types of techniques can help to determine the physical laws governing the phenomenon, giving rise to explainable artificial intelligence techniques [6].At other times, solutions are restricted to those that comply with more general physical principles, or of a higher epistemic level.Thus, inductive biases are developed to ensure (within the range of accuracy of the technique) compliance with general physical laws such as conservation [7][8][9][10][11][12], etc.
In general, physics-based machine learning aims to incorporate physical knowledge into purely data-driven strategies.The goal is to devise methods capable of learning the dynamics of physical systems from data and making accurate predictions about states not explicitly included in the training data.In the framework of physicsinformed methods, prior knowledge about the system is taken into account to improve the accuracy and interpretability of the results.Of particular interest are the so-called structure-preserving methods [13,14].In computational mechanics, geometric numerical integration refers to methods that respect the physics of a particular problem, in particular its geometric characteristics [15,16].
In particular, we are interested in learning strategies that do not take into account the exact equations that model the problem at hand, simply because it is assumed that they are not known.We focus on general formulations that only impose the fulfillment of universal physical laws, such as the laws of thermodynamics [17][18][19][20].Thus, the data serve to reveal the particular formulation of the problem under experimental conditions.Then, for new, unseen situations, the structure is predicted by appropriate interpolation.
For closed isolated systems, the General Equation for the Non-Equilibrium Reversible-Irreversible Coupling, GENERIC, formulation [21,22] defines a general structure for the evolution of the system, while ensuring the first and second laws of thermodynamics are satisfied.Numerical learning strategies based on the GENERIC formalism have successfully been developed in previous works, while numerical integrators with guaranteed stability properties also exist [23].In them, prediction is performed by appropriate interpolation on the manifold of solutions [17,24,25] or by training neural networks [26][27][28].Alternative formalisms with thermodynamical considerations, such as the generalized Onsager formalism [29,30], also offer a general structure for the evolution of reversible and irreversible processes [12,31,32].
Recently, a machine learning strategy for interacting, dissipative open systems was proposed by Hernández et al [33].It employs, during the training period, an inductive bias that generalizes port-Hamiltonian structures [34][35][36] to dissipative systems.This strategy essentially develops a port-Hamiltonian formulation, in which GENERIC is extended to open systems that communicate and exchange energy through ports.
In this paper, we present a learning strategy which applies piece-wise linear regression to the terms in the port-metriplectic formulation in [33] for interacting open systems governed by different physics.The resulting method ensures the fulfillment of the principles of thermodynamics (conservation of energy at a global scale, but exchange among sub-systems; non-negative entropy production) for the predicted states of the coupled system.The strategy is specially suitable in cases with a low amount of available data.The robustness of the method is tested for pseudoexperimental (synthetic) and experimental data.
The remainder of the paper is structured as follows.First, we briefly overview the generalized Onsager and GENERIC formalisms modeling the dynamics of closed systems in Section 2, and the port-metriplectic formulation for open systems in Section 3. Section 4 is devoted to the description of the proposed learning algorithm for the dynamics of open systems.In Section 5, we assess the performance of the strategy for three different examples.The discussion in Section 6 closes the paper.

Dynamics of closed systems
The state of a closed (isolated) system at time t is assumed to be fully described by the value of a set of variables z = z(t).We can learn the evolution of the system by identifying the structure of the dynamical problem where function f is unknown in general [37,38].
In the data-driven approach, the closed form of f is not sought.Instead, we look for an approximation of f that enables a robust and efficient integration in time with sufficient accuracy.Given some observations for the evolution of z at discrete time instants, function f is readily approximated by applying regression techniques, such as classical linear regression [17,25,39], support vector machines [40] or neural networks [26,33], among others.This classical approach to the problem presents, however, severe limitations.Very often, small perturbations in the input data produces large, physically meaningless outputs.
To void these issues, in the last years we have seen a growing interest in the development of physics-based, structure-preserving and related techniques.One particularly appealing approach to this problem assumes a particular form of f , depending on known properties of the system at hand.For instance, if the system is conservative, the Hamiltonian formalism allows us to assume that f takes the form where the Hamiltonian H is the total energy of the system, E, and L is the Poisson, skew-symmetric matrix.The problem is then reduced to determine the precise forms of L and E applying regression to data.While this type of techniques have attracted a lot of interest, purely conservative (reversible) systems are scarce in nature, where the norm is the presence of dissipation.
In what follows we briefly review two alternative formulations for the time evolution of closed systems under reversible and irreversible conditions: the generalized Onsager principle, which is a single-generator formalism, and the GENERIC framework, which is a double-generator formalism.These two approaches incorporate dissipative phenomena into the Hamiltonian formulation (2) and lead to equivalent dynamics.See [32,41,42] for a detailed discussion in the relation between single and double operator formalisms.

Generalized Onsager formalism
The generalized Onsager formalism [29][30][31] is a thermodynamically consistent formulation for the evolution of a non-equilibrium system.This is a single-generator formulation, which means that both the reversible and dissipative contributions to the evolution use the same generator, F. The generator F is a potential function with a thermodynamic interpretation such as free energy or negated entropy.The evolution of z is then modeled as where L is a skew-symmetric matrix and M is symmetric and positive semi-definite.Matrix L models the conservative part of the system, while M models the dissipative contribution.

GENERIC formalism
The General Equation for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) formalism [21,22] offers a general description for the evolution of a system based on two operators, accounting for conservative and dissipative phenomena, respectively.Within this approach, the free energy is assumed to take the form F = E+S.
With this additive decomposition, we have with L the Poisson skew-symetric matrix, modeling the conservative part of the evolution, and M a symmetric and semi-positive definite matrix, responsible for the dissipation of the system.E stands for the energy of the system and S, for its entropy.For Equation ( 4) to hold, it is necessary that the so-called degeneracy conditions hold, namely, and These conditions ensure the fulfillment of the first and second laws of thermodynamics, this is, conservation of energy and non-decrease of entropy.From ( 5) and ( 6), it is straightforward to prove given the skew-symmetry of L, and given the positive semi-definiteness of M .Thus, a regression of (1) based on the GENERIC formulation guarantees satisfying the laws of thermodynamics by construction.This approach is often referred to as metriplectic, since it combines a symplectic contribution to enforce energy conservation and a metric one to enforce entropy production [43,44].

Port-metriplectic formalism for open systems
The formalisms in Section 2 can be extended to describe the dynamics of open systems, in which the evolution of the state variables z is affected by external forces.In this framework, we use the term bulk contribution to refer to the conservative and dissipative phenomena that are inherent to the system, and the term port contribution to refer to the external forces.
In this section, we review the port-Hamiltonian extension of the GENERIC formulation to interacting and open systems proposed by Hernández et al [33], called port-metriplectic formalism.The extension of the generalized Onsager formalism to open systems is analogous [31].
The port-metriplectic formulation for the evolution of state variables z reads where we use a GENERIC description of the bulk dynamics and F accounts for the presence of an external force.In this case, satisfying the degeneracy conditions is not required: conservation of energy and increase of entropy cannot be guaranteed due to the presence of an external port.
In the work by Hernández et al [33], the port contribution F is also decomposed following the GENERIC structure.That is, with L skew-symmetric and M symmetric and positive semi-definite.E and S stand for the contributions to the energy and entropy of the system which are due to the external force, respectively.
Alternatively, we can follow the one-single operator formalism (3) to decompose F , with L skew-symmetric, M symmetric and positive semi-definite and F the port addition to the bulk free energy F. The two descriptions of the port contribution F lead to the same dynamics.

Learning algorithm
For a parametric open system, we denote by D = {z(t n ; µ j ); n = 1, . . ., n T , j = 1, . . ., n p } the dataset with observations of its state variables z at discrete time instants t 1 , . . ., t n T , for a representative parametric sampling, µ 1 , . . ., µ np .In this section, we present a strategy to approximate z(t n ; µ ⋆ ), n = 1, . . ., n T , for a parameter µ ⋆ which is not included in the original sampling.The strategy is based on interpolating the terms in the port-metriplectic formulation (9) of the problem.We distinguish two phases in the process: identification and interpolation.The identification phase consists in the computation of all the elements in (9) (matrices and gradients), at every time step and for all the sampled parameters in the data set.Then, in the interpolation phase, the precomputed elements are properly combined to approximate the port-metriplectic structure for the new parameter µ ⋆ .The details of each phase are discussed next.

Identification (offline computation)
The first step is computing the value of the matrices and gradients appearing in the portmetriplectic structure (9) of the problem for the parameters µ j , j = 1, . . ., n p , in the data set D.
Note that this is an offline computation, to be performed once at the beginning of the process.At time step n, the discretization of Eq. ( 9) leads to where we use the forward Euler approximation of the derivative.Although more accurate integration schemes can be employed to minimize the loss of accuracy during the integration, the study of the effects of these schemes in the resulting prediction is not the objective of this work.For a detailed analysis on the topic, the interested reader is referred to [23].
Assuming that the energy and the entropy depend quadratically on the state variables, we approximate where A n and B n are assumed to be diagonal matrices.Of course, more general assumptions can be made, as in [17], but for the examples to follow, this mild simplification has rendered excellent results.
In the following, we take a GENERIC decomposition of the port term F n as in (10), that is, we express with A n and B n diagonal matrices.The formulation for the single-generator decomposition of the port ( 11) is analogous.
The discrete values of the terms in the equation are computed by solving the minimization problem where we have omitted the dependence on µ j to simplify the notation.Very often, the Poissson and dissipation matrices have a known structure that can be fixed beforehand, and only the gradients (or equivalently, matrices A and B) remain as unknowns.

Staggered approach
In cases in which data of the equivalent closed system (with no external force acting) are available, the minimization problem in (15) can be solved by means of a staggered approach, in order to facilitate its convergence.In particular, the GENERIC formulation (4) of the closed system is used to determine the bulk contribution of the portmetriplectic open-system formulation.Let C = {z c (t n ; µ j ); n = 1, . . ., n T , j = 1, . . ., n p } denote the dataset containing the observations for the closed system.In the staggered approach, we take the next two steps: 1. Computation of the terms from the bulk contribution, using the closed-system data in C.
The evolution of the closed system is captured by the GENERIC formulation (4) of the problem.The time discretization of (4) leads to the constrained minimization problem subject to 2. Computation of the port contribution, using the open-system data in D. The discrete terms in the bulk contribution are replaced by the values obtained in (16).Thus, the remaining terms in F n+1 are obtained by solving Note that we are extrapolating the bulk opensystem gradients DE n and DS n from the closedsystem gradients, by using the matrices A n and B n that are obtained from the data in C.

Interpolation (online computation)
For a new parameter µ ⋆ , we approximate the time evolution of the corresponding state variables, z(t n ; µ ⋆ ), using the discrete port-metriplectic formulation (12).That is, given the initial state variables z 1 , for n = 1, . . ., n T − 1, with the matrices in the expression to be determined.
Here, we propose to approximate the unknown matrices in the formulation by linear interpolation on the manifold of GENERIC terms.Let N denote the set of indices of the neighboring parametric solutions in D, enclosing the value of µ ⋆ in the parametric space.Recall that the portmetriplectic terms corresponding to time evolution for µ i , i ∈ N , are precomputed and stored in the identification phase.Then, at time step n, we approximate for some weights w i , i ∈ N .The approximation for all unknown matrices in (18) in terms of precomputed matrices for neighboring solutions is analogous.
The weights w i are here computed from a linear interpolation in the parametric space of µ, with w i ≥ 0 for all i ∈ N and i∈N w i = 1.That is, the regression is performed globally for all parameters.

Numerical results
In what follows, we show the performance of the methodology proposed in Section 4 in three different examples involving different physics.Recall that we distinguish two phases: identification and interpolation.In the identification phase we approximate the discrete port-metriplectic structure (12) for a given dataset.In the interpolation phase, we use the learned numerical integrator to predict the evolution of the system in unseen scenarios.

Damped harmonic oscillator
We consider a one-dimensional harmonic oscillator in the presence of friction [45].The system can be interpreted as an open system, in which a perfect harmonic oscillator (bulk) is damped because of the action of an external force (port).
We take the independent state variables z = (q, p, S), with q the position of the particle, p its momentum and S the entropy of a homogeneous medium causing the friction of the particle.The motion is described by where m is the mass of the particle, T is the temperature, k is the spring constant and γ is the damping coefficient.The initial condition is of the form (q 0 , 0, 0).We set m = 1 kg, T = 25 K and k = 2250 N/m.The evolution of the state variables can be modelled by a GENERIC structure (4), with known Poisson and dissipative matrices [45].In particular, Here, we aim to model the system by the portmetriplectic formulation (9), through its discretized version (12).The Poisson and dissipative matrices from both the bulk and port contributions are predefined as where we take the Poisson matrix of the GENERIC formulation, and for the dissipative matrix, we only indicate the zero coefficients and let the gradient matrices B and B acquire the constants.
We account for two different datasets, obtained by integrating the equations for (i) 10 different values of the damping coefficient γ, uniformly distributed in [0.2, 2] Ns/m, and fixed initial position q 0 = −0.075m, and (ii) 11 different values of the initial position q 0 ∈ [−0.15, −0.075] m and fixed γ = 1 Ns/m.The snapshots in the datasets are depicted in Figure 1.
The time discretization is uniform in the interval [0, 10] s, with time increment ∆t = 10 −3 s.The evolution of the corresponding closed system (undamped harmonic oscillator) is also integrated for all cases, with γ = 0.The terms in the port-metriplectic formulation are then obtained by following the staggered approach described in Section 4.1.
To test the learning algorithm, we interpolate the intermediate elements of the datasets using the GENERIC terms arising from two neighboring snapshots.Since the parametric discretizations are uniform, we have w 1 = w 2 = 0.5 in the linear combinations in (19).The obtained errors are summarized in Figure 2.For the first dataset, with parametrized damping coefficient, we obtain relative L2 errors around 10 −2 .For the second one, with parametrized initial position, errors are of the order of 10 −3 .

Vertical sloshing tank
In this example, we use experimental data from the SLOWD database [46,47].The experimental setup is a liquid-filled tank, which is attached to a spring-damper system.Figure 3 shows an scheme of the setting.This experiment was thought to characterize the effects of the fuel sloshing in airplane tanks.This fuel acts as a dampener for the wing vibration.
The structure is deflected for the tank to start the motion from an initial position q 0 .Then, the tank is released and starts oscillating in the vertical direction.The liquid that fills the tank acts as an external force that interacts with the dynamics.
The available data contains measurements of the vertical position and acceleration of the tank and the load cell force for several configurations.As a post-process, we are able to obtain additional variables such as the momentum, the entropy, and the sloshing force.For the specifications on the setup and the measuring system we refer to Martínez-Carrascal and González-Gutiérrez [46].
The chosen state variables to describe the system are z = (q, p, S), with q the vertical position of the tank, p its momentum and S the entropy.
We test the proposed learning methodology for two different datasets.In the first set, we account for 9 different filling levels of liquid, from 10% to 90% of the capacity of the tank, starting the motion from the position q 0 = −0.075m.The time interval is [0, 6.5] s and it is discretized with time increment ∆t = 9.8 • 10 −4 s.
In the second set, we consider the observations for 10 different initial positions q 0 distributed in the interval [−0.055, −0.011] m, and a fixed filling level of 50%.The time interval is [0, 4.5] s and it is discretized with time increment ∆t = 9.4 • 10 −4 s.The position evolution for the experiments included in the datasets are shown in Figure 4.
For all considered configurations, measurements for the equivalent experiments with the dry tank are also available.This enables the use of the staggered scheme in the identification phase of our algorithm.Due to the oscillating behavior of the observations, the Poisson and dissipative matrices in the port-metriplectic formulation are predefined as in Eq. ( 22).
Fig. 2 Damped harmonic oscillator.Boxplots for the relative L2 errors in the interpolation of all the snapshots in the two datasets with (i) parametrized damping coefficient γ, and (ii) parametrized initial position q 0 .The state variables are position q, momentum p and entropy S.
Fig. 3 Vertical sloshing tank.Experimental setup for the data in the SLOWD database [46,47].The setup consists of a rigid tank, attached to the walls by two spring-damper systems.One of them is also equipped with a load cell, not represented for simplicity.
For previously unseen situations (filling level or initial position as parameters), interpolation is performed using the GENERIC constituents coming from two neighboring snapshots.Figures 5 and 6 show the state variables position and entropy for previously unseen situations using the two datasets.Time integration is performed in the whole time intervals, and is zoomed at [1.6, 1.9] s for illustration purposes.The learned solutions are able to capture the oscillating behavior of the system and, as expected, are located somewhere between the neighboring solutions.This motivates the use of enclosing solutions as neighbors.Figure 7 summarizes the L2 relative errors for interpolation of all available experiments.In both cases, the obtained errors are slightly larger than the errors for the synthetic damped harmonic oscillator example in Section 5.1.However, the mean of the errors is below 10% even in the presence of noise in experimental data.
The strategy performs robustly despite the presence of experimental noise in the data.It is worth mentioning that a naive approach based on pure interpolation of the results by taking the amount of filling of the tank as the governing variable does not lead to any meaningful result.

Fluid-structure interaction in an oscillating tank
In this case we analyze data coming from a fluidstructure interaction problem taken from [48].In it, see Fig. 8, a cylindrical tank is attached to a cantilevered beam.Both the solid and fluid dynamics in this experiment share similar characteristics, leading to strong couplings.In the experimental setting, two piezoelectric actuators are attached to the beam, near the encastre.These cause the beam to oscillate in torsion.The motion of the set is captured with the help of two accelerometers at the free edge of the beam.The movement of the tank in the horizontal direction causes the beam to bend along its weak axis of inertia, while the weight of the tank and the liquid make it bend in the vertical direction.Additionally, oscillations of the tank cause torsion in the beam.On top of all these physical phenomena (two bending directions, torsion and sloshing;   the tank is assumed to be perfectly rigid), sloshing of the fluid dampers the dynamics of the tank.The tank is 0.5 m wide, with an internal cavity of 0.470 m and an internal diameter of 0.105 m.It is made of a material with mass density 1180 kg/m 3 and filled with water.The beam is assumed to be made of aluminum, with Young's modulus 75 GPa, Poisson's ratio 0.33 and mass density 2970 kg/m 3 .Its dimensions are 1.36 × 0.15 × 0.005 m 3 .
To obtain synthetic data, the same method developed in [48] is employed.It uses a port-Hamiltonian approach and an appropriate discretization for each physics.The interested reader is referred to this article for more details about the particular implementation.In essence, the beam is assumed to follow a classic Euler-Bernoulli-Navier model for the bending phenomenon, the fluid is assumed to follow the shallow-water equations and, finally, torsion is assumed to be of Saint Venant type.The tank moves as a rigid body.Numerical results obtained under these assumptions are assumed to be the ground truth for our method.
Our vector of (synthetic) measurements is composed by z = [x B , x T , x F , x RB ] ⊤ , where The actuation on the beam causes the tank to rotate (depicted in dashed line) but also to oscillate in the horizontal direction.Its weight also causes the beam to bend in the vertical direction (not represented for simplicity).The movement in the tank causes the fluid to slosh.
• x RB ∈ R 6 represents the degrees of freedom of the rigid-body motion of the tank (three displacements, three rotations of the center of mass).
In the simulations taken as ground truth, N B = N T = N F = 10.
The system is analyzed for different degrees of filling of the tank, following a staggered approach, as introduced in Section 4.2.Two different data sets are considered, one with 40, 50, . . ., 80% of filling, and a more detailed one, with 40, 45, 50, . . ., 80% of filling.Following the staggered approach, we first determine the time evolution of each sub-system by identifying the elements of their GENERIC description, and then we identify the contribution of each port, The error in the reconstruction of a previously unseen degree of filling is evaluated as the ℓ 2 -norm error in the fluid surface.Thus, assuming that there are n nodes discretizing the fluid surface, this error is computed as where z GT refers to the height at a given nodal position in the ground truth and z app its approximated counterpart, provided by our method.
With these settings, we obtained the errors reported in Table 1 for different filling levels.The reader can notice the excellent degree of accuracy obtained for both datasets.In practice, it seems that the refined dataset does not always provide a substantial increase in accuracy.In order to ascertain if a more elaborated interpolation scheme could provide with more accurate results, we have also tested a quadratic scheme, that employs three neighboring snapshots to determine the GENERIC structure of the problem.As an example, the 61 % filling provides 0.487036 % error for the coarse dataset and 0.411062 % for the detailed dataset.
The evolution in time of the error is depicted in Fig. 9 for the 47 % of filling example.
Finally, for a qualitative evaluation of the error, some snapshots are provided in Fig. 10.

Conclusions
We developed a method able to provide with accurate estimates of the dynamics of coupled, multi-physics and parametric systems from data.The method is based first on the regression from experimental or synthetic data of the terms of an assumed port-metriplectic structure for the problem at hand.This ensures that the learnt evolution of the system will satisfy by construction the first and second principles of thermodynamics, even for dissipative, open systems.Then, in a second step, for previously unseen situations, the method interpolates each term of the metriplectic description from neighboring parametric data.Unlike previous approaches, the just developed method does not employ neural networks, that have demonstrated to be a powerful tool, but classical constrained regression techniques.It nevertheless offers competitive results.The method is robust for synthetic as well as experimental (noisy) datasets, and is able to provide with complete rollouts of the evolution of the different systems tested until stop by dissipation.
It remains as a topic for future analysis wether the employ of existing neural network architectures could offer competitive advantages over the just presented constrained regression approach.In any case, the assumed port-metriplectic structure of the evolution of the different considered systems has shown to be a powerful inductive bias for learning complex, multi-physics systems from data.

Fig. 6
Fig. 6 Vertical sloshing tank.(ii) Parametrized initial position.Learned solution for initial position of −0.0334 m, interpolating the structures of neighboring snapshots with initial positions −0.0281 m and 0.0384 m.Plot of the position q and entropy S at time interval [1.6, 1.9] s.

Fig. 7
Fig.7Vertical sloshing tank.Boxplots for the relative L2 errors in the interpolation of all the snapshots in the two datasets with parametrized (i) liquid filling level, and (ii) initial position.The state variables are position q, momentum p and entropy S.

Fig. 8
Fig.8Oscillating tank.Sketch of the experimental setup.The actuation on the beam causes the tank to rotate (depicted in dashed line) but also to oscillate in the horizontal direction.Its weight also causes the beam to bend in the vertical direction (not represented for simplicity).The movement in the tank causes the fluid to slosh.

Fig. 9
Fig. 9 Oscillating tank.Evolution in time of the L2-norm error of the position of the free surface of the fluid.

Table 1
Oscillating tank.Errors in L2 norm in the approximation of the liquid free surface by employing the detailed and coarse datasets, respectively.