Michael Kirsten
email: kirsten@kit.edu

Olivier Cailloux
email: olivier.cailloux@dauphine.fr

Towards automatic argumentation about voting rules

Keywords: Social choice theory, bounded model checking, argumentation theory, automated reasoning

Voting rules aggregate the preferences of a group to make decisions. As multiple reasonable voting rules exist, the axiomatic approach has been proposed to exhibit both their merits and paradoxical behaviors. It consists in characterizing a voting rule by a set of understandable properties called axioms. It is however a difficult task to characterize a voting rule by such axioms, and even when a proof exists, it may be difficult to understand why a specific voting rule fails to satisfy a given axiom, especially for untrained users. In this article, we present an automatic method which determines whether a given rule satisfies a set of axioms. When the rule does not satisfy an axiom, the automatic prover generates comprehensible evidence of the violation in the form of a counter-example. It can be used by non-expert users to comprehend the violation and may serve to argue in favor of other rules which satisfy the axiom. Our method is based on the software analysis technique bounded model checking, which enables bounded verification for properties of software programs. It translates the program together with user-annotations into a reachability problem for those profiles and outcomes which adhere to our specification. The method can be applied to arbitrary voting rules; we demonstrate it on the case of the Borda axiomatization and compare the Borda rule to both the Black and the Copeland voting rules.

Goal and motivation

Voting rules serve to select winning alternatives from a set of a priori possible alternatives, given the preferences of a set of voters. Many voting rules that appear reasonable have been proposed in the literature, however each of them shows paradoxical behavior for at least some voting situations [START_REF] Arrow | Social choice and individual values[END_REF]. Axiomatic analysis permits to distinguish the merits of different voting rules: it consists in determining axioms, properties of voting rules that seem a priori desirable for a voting rule to have, and sets of axioms which together characterize a voting rule. A set of axioms characterizes a voting rule when there is a unique voting rule that satisfies all axioms in the set.

Our first interest in this article is, given a voting rule, to detect automatically which axioms it satisfies among a given set of axioms. Although it is not possible to do this efficiently in full generality, modern SAT-and SMT-solvers [START_REF] Biere | Symbolic model checking without BDDs[END_REF] enable solving this task efficiently when limiting the number of voters and alternatives. Indeed, as stated in the "small scope hypothesis" and also seen in related work on automated reasoning of voting rules, a high proportion of violations can often already be found for a comparatively small fraction of voters and alternatives [START_REF] Jackson | Software Abstractions -Logic, Language, and Analysis[END_REF][START_REF] Brandt | Optimal bounds for the no-show paradox via SAT solving[END_REF]. With this prerequisite, the technique called software bounded model checking (SBMC) allows encoding the problem into a decidable set of logical equations, which can then be solved by a SAT-or SMT-solver. In this article we propose an automatic method that determines whether a rule satisfies an axiom within the specified bounds. We illustrate it by applying the approach to the case of Borda: we consider a set of four axioms that characterize the Borda rule, and given some rules different than the Borda rule, we find automatically which axioms from that set each rule fails to satisfy.

It is especially interesting to obtain short and readable proofs that a rule does not satisfy some axiom. This is our second objective in this article: our automatic prover outputs profiles, preferably small or satisfying some other property representing "simplicity" for a human eye, that exhibit a failure of a rule to satisfy an axiom. Such a short proof is interesting even if it is already known in the literature that the rule does not satisfy the axiom, as short proofs permit to help getting an intuitive understanding of voting rules. Furthermore, it can be used to talk about voting rules with non-expert users. Our illustration in the case of Borda permits to obtain a tool that "argues in favor" of Borda [START_REF] Cailloux | Arguing about voting rules[END_REF] (assuming the axioms characterizing Borda of our choice are considered desirable properties of any voting rule): given any rule that is not the Borda rule, the prover finds a concrete profile which illustrates that the rule fails to comply with some of those axioms.

We consider an infinite set of voters N and a fixed finite set of alternatives A . An election involves a subset of voters N ⊆ N that each express preferences over all alternatives. Those preferences are elements of L (A ), the linear orders over the alternatives: connected, transitive, asymmetric binary relations over A . Given a voter i ∈ N , we write as ≻ i ∈ L (A ) the preference of voter i. Given a set of voters N, a profile (≻ i ) i∈N = R is an association of each voter to her preference. The set of all possible profiles is N⊆N L (A ) N . A voting rule f : N⊆N L (A ) N → P * (A ), where P * (A ) represents the non-empty subsets of A , maps each possible profile (≻ i ) i∈N ∈ N⊆N L (A ) N to a non-empty set of winning alternatives / 0 ̸ = W ⊆ A (sets of winners are used to account for tied winners). An axiom is here simply defined as a property of a voting rule. One kind of axioms of interest (for their simplicity) are so-called functional properties [START_REF] Beckert | Automated verification for functional and relational properties of voting rules[END_REF]. A functional property is one that can be represented as a relation S ⊆ N⊆N L (A ) N × P * (A ): given any

R ∈ N⊆N L (A ) N , S(R) = { A | (R, A) ∈ S }
contains all the admissible sets of winners given that profile. When S(R) = P * (A ), S expresses no constraint about that profile. A voting rule f satisfies the functional property S iff f picks sets of winners among the ones considered possible by S, namely, iff f ⊆ S (where we view both f and S as relations). Whenever a rule f fails to satisfy a functional property, it can be proved in a very short and understandable way, by exhibiting a profile R on which f picks a wrong set of winners, f (R) / ∈ S(R). The property DOM is an example of a functional property.

DOM This property forbids Pareto-dominated alternatives to win. An alternative a ∈ A is Pareto-dominated in a profile

(≻ i ) i∈N iff some alternative a ′ ∈ A is unan- imously preferred to a in (≻ i ) i∈N : ∀i ∈ N, a ′ ≻ i a.
The property mandates that the rule selects winning alternatives among U R , denoting the alternatives that are not Pareto-dominated in R:

∀R ∈ N⊆N L (A ) N , S DOM (R) = P * (U R ).
Another, more general, kind of axioms we are interested in, are k-relational properties, with k ∈ N [START_REF] Beckert | Automated verification for functional and relational properties of voting rules[END_REF]. A k-relational property S ⊆ ( N⊆N L (A ) N × P * (A )) k contains all the admissible ways of associating sets of winners to k profiles. A voting rule f satisfies a k-relational property S iff f k ⊆ S, viewing again f as a relation. For example, for

k = 2, if S contains (R 1 ,W 1 , R 2 ,W 3 ) and (R 1 ,W 1 , R 2 ,W 4 ), given a voting rule f such that f (R 1 ) = W 1 , S declares as admissible that f (R 2 ) be equal to W 3 or W 4 .
The property REINF exemplifies a 3-relational property.

REINF Reinforcement requires that elections uniting disjoint groups of voters elect the candidates chosen by both groups, if some such candidates exist: for each

(≻ i ) i∈N 1 , (≻ i ) i∈N 2 with N 1 ∩ N 2 = / 0 and A ̸ = / 0, defin- ing A as f ((≻ i ) i∈N 1 ) ∩ f ((≻ i ) i∈N 2 ), reinforcement im- poses that f ((≻ i ) i∈N 1 ∪N 2 ) = A.

Checking properties automatically

Let us start by describing the general technique of checking for software properties, before we explain how we apply it to verify whether voting rules satisfy properties. In a nutshell, we translate the proving task to a reachability problem, which the SBMC tool encodes into an instance of the SAT problem.

Checking software properties

We assume we are given an algorithm A under the form of an imperative program (for example, written in the C language), that uses some parameter values taken among a set of possible values I. An entry i ∈ I is a list of values, one value for each such parameter: it gives a value to everything that a run of A depends on, such as its input variables, or anything that is considered non-deterministic from the point of view of A. For this reason, those parameters are qualified as "non-deterministic", to distinguish them from normal parameters used in a programming language to pass information around. (By contrast, some values can be "derived", thus, computed in A from the non-deterministic parameter values, or declared as constants in A, and both values of non-deterministic parameters or derived values can then be used as normal parameters in the program.) We are also given a software property to be checked about A, in the form C ant ⇒ C cons , where ant and cons stand for antecedent and consequence respectively. Both C ant and C cons are sets of boolean statements. A boolean statement is a statement of A that evaluates to a boolean value, for example, a statement checking that some computed intermediate value is odd. An entry i is said to satisfy a set of boolean statements iff all boolean statements in the set evaluate to true during the execution of A using the non-deterministic parameter values i, and is said to fail the set of boolean statements otherwise. The property C ant ⇒ C cons requires that for all possible entries i ∈ I, if i satisfies C ant , then i satisfies C cons . As an example, assume A computes, given i, two intermediate integer values v 1 and v 2 , and then returns a third value v 3 . The property to be checked could be: if v 1 is negative, then v 2 is positive and v 3 is odd. A solver that is asked to check a software property C ant ⇒ C cons thus exhaustively searches for an entry i that satisfies C ant but fails C cons . The property is valid iff it is impossible to find such an entry. SBMC is a fully-automatic static program analysis technique used to verify whether such a software property is valid, given an algorithm and a property to be checked. It covers all possible inputs within a specified bound. It is static in the sense that programs are analyzed without executing them on concrete values. Instead, programs are symbolically executed and exhaustively checked for errors up to a certain bound, restricting the number of loop iterations to limit runs through the program to a bounded length. This is done by unrolling the control flow graph of the program and translating it into a formula in a decidable logic that is satisfiable if and only if a program run exists which satisfies C ant and fails C cons . The variables in the formula are the non-deterministic parameters of A, and their possible values are taken from I. This reduces the problem to a decidable satisfiability problem. Modern SAT or SMT-solving technology can then be used to verify whether such a program run exists, in which case an erroneous input has been found, and the run is presented to the user. If the solver cannot find such a program run, it may be either because the property is valid, or because it is invalid only for some run which exceeds the bound. In some cases, SBMC is able to infer statically which bound is sufficient to bring a definitive conclusion.

Checking voting rule properties

In our case, we are given a voting rule f implemented as a C function that returns winners given a profile, and a list of axioms. When the axiom is a functional property S, in order to check whether f satisfies S, we want to use SBMC to search for an input profile R such that f (R) / ∈ S(R). To do this, we implement the requirements of S as a C function with annotations that indicate to the SBMC tool a set of boolean statements, constituting C ant , that evaluate to true iff f (R) / ∈ S(R). Other annotations serve to indicate nondeterministic parameters, including the input profile, and other non-deterministic parameters depending on what is necessary for our implementation to check whether a profile satisfies S. We use as C cons a single boolean statement that is always false. This makes the solver search for an entry R (and possibly other parameter values) such that f (R) / ∈ S(R). For example, when searching for a violation of Paretodominance (DOM as defined in Section 2), it suffices to implement a C function F containing boolean statements that evaluate to true whenever some alternative a, a nondeterministic parameter of F, is Pareto-dominated by another alternative a ′ , another non-deterministic parameter of F, but wins. This way, the solver will search for an entry (R, a, a ′ ) corresponding to an input profile that contains a Pareto-dominated alternative that wins, thus, an input that violates S, and will report such an input, if found. This approach can be extended to cope for k-relational properties, where k > 1. We do not detail this as this is similar in spirit, only requiring to duplicate the input parameters so that the solver searches for multiple profiles as entries, and writing the function checking for S so that it cares about the relationship between multiple sets of profiles and winners. In order to apply SBMC, we need to set bounds for the number of loop iterations. We have implemented all tested axioms and voting rules (and more) in such a way that the number of loop iterations depend only on the number of voters and alternatives. In the experiments, we set those as constants. The SBMC tool can thus automatically infer the right bounds in a syntactic pre-processing step. For more complex cases where it would not be possible to syntactically infer the bounds, e.g., when while-loops are used or the number of iterations involves some non-deterministic parameter, the bound could be specified manually. Regarding computation time, finding an input that violates S, if such an input exists, scales much better than proving that no such input exists, as we are allowed to stop as soon as such an input is found. Moreover, our method generally scales better for smaller numbers of loop iterations, simply based on the smaller formula for which a model is searched. As we are generally looking for "small" profiles, this is most often the case. For very complex (i.e., having complex loop structures with a large number of iterations) voting rules and axioms, or for properties which are not violated, it may happen that the solver takes more computation time or resources than the user is willing to invest. Our method is generally well-suited for finding counter-examples to voting rule properties as most violations are exhibited for relatively small bounds. We substantiate these claims in the following section, where we exemplarily use our method to argue in favor of the Borda rule. This approach can be generalized to search for profiles that illustrate the satisfaction of simultaneous axioms, or profiles that satisfy some axioms and not some others, etc. We have written a general program to search for such situations that cover about fifteen classical axioms 1 . For our experiments, we use CBMC 5.8 [START_REF] Clarke | A tool for checking ANSI-C programs[END_REF], an implementation of the SBMC approach for the C language, with the built-in solver based on the SAT-solver MiniSat 2.2.0 [START_REF] Eén | An extensible SAT-solver[END_REF]. All experiments are performed on an Intel(R) Core(TM) i5-6500 CPU at 3.20 GHz with 4 cores and 16 GB of RAM.

A simple example

Listing 1 illustrates how we specify DOM, the "Paretodominance" functional property. The function is annotated with specifications in the form of assumptions (expressed by the function __CPROVER_assume), that represent the statements composing C ant , as well as operations to model non-deterministic choice (indicated by the prefix nondet_) that represent what we called non-deterministic parameters in Section 3.1. These non-deterministic choice operators can be used anywhere inside the analyzed program and are also translated to the formula which is passed to the solver. The solver then searches for instantiations of these variables which lead to a violation of C cons . In the displayed function, the constants N and M denote the numbers of voters and alternatives respectively. The voters and alternatives are represented by integers from 0 to N -1 and from 0 to M -1. The function accepts as input a profile prof, modeled as a two-dimensional array:

prof[i][r]
is the alternative that voter i associates to rank r, where 1 The runnable source code can be found at https://github.com/mi-ki/ voting-rule-argumentation. 0 is the best, and M -1 the worst rank. The function also accepts as input a set of winning alternatives, modeled as an array win with one entry per alternative: win[a] equals 1 if the alternative a belongs to the set of winning alternatives, 0 otherwise. The non-deterministic parameters are prof (whose declaration is not shown here) and the two alternatives bad and good (Line 2). Lines 6 to 14 initialize the array prefergtob so that prefergtob[i] holds the value 1 if voter i prefers good to bad, 0 otherwise. The function thus indicates to the solver that it must find two alternatives bad and good, thus integers in the suitable range (Lines 3 and 4) and different from each other (Line 5), such that everybody prefers good to bad (Lines 15 and 16), and yet the alternative bad is a winner of the election (Line 17). Any run which satisfies these specified statements is a valid counter-example which would prove that the profile and winners given as input violate Pareto-dominance. We pursue the example with the setup given in Listing 2. Therein, we initialize a profile with symbolic nondeterministic values (Line 7) and restrict their ranges such that they are valid alternatives (Line 8). Furthermore, we use a helper array used with one entry per alternative, which is used to ensure that every ballot holds every alternative only once (Lines 5, 9 and 11). We assume a voting rule is given as a function f (Line 14) which gets a profile prof and the number of voters N as parameters, and returns the set of winning alternatives as output. For the experiments within this article, the implementations of f will follow directly from their definitions. After calling the test methods for the properties to be checked (in this case only Pareto-dominance), we set a boolean statement that is always false (Line 16) as content of C cons (recognized by the solver by the keyword assert), which indicate to the solver that any program run that reaches this point is a counter-example of interest to us.

Definitions for the experiments

We propose here to consider the axioms that characterize the Borda voting rule (Def. 1), as defined in a variant of Young's axiomatization [START_REF] Cailloux | Arguing about voting rules[END_REF]. We also define a few axioms that are not satisfied by the Borda rule, and finally define two supplementary rules. All these concepts will be used in the experiments. (Readers knowledgeable in social choice theory can skim this section.) Definition 1 (Borda). The Borda rule, given a profile (≻ i ) i∈N , associates to each alternative a and voter i the score s(a, i) equal to the number of alternatives that a beats in ≻ i , and associates to each alternative a the score s ′ (a) = i∈N s(a, i). The winners are the alternatives that have the maximal score: f Borda ((≻ i ) i∈N ) = arg max a∈A s ′ (a).

An axiomatization of the Borda rule

We first require a few definitions. A profile (≻ i ) i∈N is elementary iff it has exactly two voters, and if the set of alternatives can be partitioned into disjoint subsets A top , A bottom ⊆ A with A top ∪ A bottom = A and A top ̸ = / 0 such that both voters have all alternatives in A top preferred to all alternatives in A bottom , and the voters have inverse preferences over A top , and inverse preferences over A bottom . Thus, denoting the voters by 1 and 2, and denoting by ≻ | A the restriction of ≻ to A, ∀A ∈ {A top , A bottom } : We say that a profile (≻ i ) i∈N is cyclic iff it has exactly m voters and m different linear orders, and some cycle corresponds to all its linear orders (equivalently, the linear orders in (≻ i ) i∈N are all those that correspond to a given cycle).

≻ 1 | A = (≻ 2 | A ) -1 . Let T (R)

△

Below is the axiomatization that we use for the Borda rule, composed of the axioms ELEM, CYCL, CANC and REINF. It is very similar but not identical to the axiomatization given by [START_REF] Young | An axiomatization of borda's rule[END_REF]. The proof that these four axioms characterize the Borda rule is given in [START_REF] Cailloux | Arguing about voting rules[END_REF].

ELEM This axiom mandates that the rule, when given any elementary profile, selects its top alternatives:

S ELEM (R) = {T (R)} if R is an elementary profile and S ELEM (R) = P * (A ) otherwise.
CYCL This axiom requires the rule to select all alternatives as tied winners when given any cyclic profile: S CYCL (R) = {A } if R is a cyclic profile and S CYCL (R) = P * (A ) otherwise.

CANC The cancellation axiom constrains the set of winners to be A when all pairs of alternatives (a, b) are such that a is preferred to b for as many voters as b is to a in (≻ i ) i∈N .

REINF As defined above.

Two axioms not satisfied by Borda

Borda notoriously fails (when |A | ≥ 3) to satisfy the following two functional properties COND and MAJ, as well as the derived property WMAJ. Given a profile R = (≻ i ) i∈N , we say that an alternative a obtains a strict majority against a ′ , denoted by a M R a ′ , iff more than half of the voters prefer a to a ′ in WMAJ This property requires that, whenever some alternative is placed first by more than half of the voters in (≻ i ) i∈N , it becomes a (not necessarily unique) winner.

(≻ i ) i∈N : a M R a ′ ⇔ |{ i | a ≻ i a ′ }| > |{ i | a ′ ≻ i a }|.

Two Condorcet compatible voting rules

To end this section, we define two famous voting rules that satisfy the Condorcet property.

Definition 2 (Black). The Black [1958, p. 66] rule selects the Condorcet winner if there is one, otherwise, the Borda winners.

Given a profile R = (≻ i ) i∈N , let M R (a) denote the set of alternatives against which a obtains a strict majority, and M R -1 (a) the set of alternatives that obtain a strict majority against a.

Definition 3 (Copeland). The [START_REF] Copeland | seminar on applications of mathematics to social sciences[END_REF] rule (actually a close variant of a rule proposed by Ramon Llull in the 13th century [START_REF] Colomer | Ramon Llull: from 'Ars electionis' to social choice theory[END_REF]), given a profile R, gives to each alternative the score s(a

) = |M R (a)| -|M R
-1 (a)|, and lets the alternatives with maximal score win.

Experiments

In this section we first experiment with the Borda rule itself. In Section 5.1, we look at whether the solver is able to find out that the Borda rule satisfies Pareto-dominance, depending on the bounds we set on the numbers of alternatives and voters. We then check how long it takes to find an example that illustrates Pareto-dominance (one where some alternative is dominated and indeed not included among the winners). In Section 5.2, we search for counter-examples that illustrate that the Borda rule fails to satisfy the properties defined in Section 4.2. Second, we experiment with the Borda axiomatization. We suppose we are given a rule f , different from Borda (Def. 1), as a C function: in our experiments we consider the Black rule and the Copeland rule. In this simple experimental setup, we thus know that f fails to satisfy at least one of the Borda axioms from Section 4.1. For each of these rules, we illustrate that we can find out automatically which axioms f fails to satisfy, and output a short proof of this, easy for a human to inspect. We also analyze in which situations we can prove that axioms are satisfied by f , for those that f satisfies.

Borda and Pareto-dominance

In Figures 1 and2, we illustrate the sizes of situations, i.e., numbers of voters (on the x-axis) and alternatives (each having a different plot), that we are able to analyze. In this example, as well as in the following experiments, we only consider run-times below 30 minutes to be reasonable, and stop computations which require more time than 30 minutes. Run-times are given in seconds (on the y-axis) for up to 30 voters and 5 alternatives. Figure 1 shows the run-times for the verification that Pareto-dominance (Listing 1) holds for the Borda rule (the line for 2 alternatives is almost superposed to the x-axis). Experiments for 5 alternatives and more than 26 voters took more than 30 minutes and are thus not plotted. Figure 2 shows the run-times for finding an example situation that illustrates that Borda satisfies DOM. We used the code from Listing 1 with the only difference that we negated the statement in Line 17. The function now indicates to the solver that it must find two different alternatives bad and good, such that everybody prefers good to bad, and the alternative bad is not a winner of the election. In the first experiment, the solver proves that no profiles satisfy the provided conditions. We can see that the verifi- cation is feasible for small sizes of profiles within 15 minutes. Expectedly, in situations where some profile exists that satisfies the provided conditions, the run-times are significantly smaller. They stay well below 10 seconds for this example. Both experiments indicate that our method yields reasonable run-times for at least up to five alternatives and 25 voters.

Counter-examples to Borda

We illustrate here our approach by finding counterexamples which show that the Borda rule does not satisfy the properties defined in Section 4.2 (those properties are satisfied by both the Black and the Copeland voting rule, however). We start with the stronger Condorcet property COND. The smallest example proving that Borda fails COND can be found in less than one second for three voters and three alternatives:

Automatic comparison of Borda with other voting rules

We here want to compare the Black rule and the Copeland rule with the axiomatization of the Borda rule shown before. We first encode all four axioms in a way similar to the dominance property shown in Listing 1. The axioms ELEM, CYCL and CANC are functional properties and can thus be encoded in the very same manner. As the axiom REINF is a 3-relational property, it needs further statements relating the profiles. The three profiles are initialized with symbolic non-deterministic values as seen in Listing 2. We consider the bound N as the number of voters of the joined profile prof, and fix the sizes of the two subprofiles prof1 and prof2 using a non-deterministic value s to define prof1's size as s, and the size of prof2 as N -s. We furthermore fix the profiles prof1 and prof2 with respect to prof using an array with non-deterministic values to model a mapping between the joined profile and its two sub-profiles.

Having encoded all four axioms as either functional or krelational properties in C functions, we can now compare the Black rule and the Copeland rule with this axiomatization. We hereby rely on the "small scope hypothesis" and focus on small profiles with up to five alternatives and eight voters.

Our method successfully verifies that both the Black rule and the Copeland rule satisfy the axiom ELEM in less than 17 seconds (for "small profiles" as defined above). As this axiom must only be verified for two voters, only the number of alternatives must be bounded.

For the axiom CYCL, the verification for both voting rules can be done within 10 minutes (also addressing only "small profiles" as defined above). Here, the numbers of voters are fully determined by the numbers of alternatives (CYCL only addresses profiles with an equal number of alternatives and voters).

When checking the axiom CANC, for both rules, verification is achieved within a few minutes for the cases of eight voters or less when the number of alternatives is bounded to three; five voters or less when the number of alternatives is bounded to four; and three voters or less when the number of alternatives is bounded to five. As run-times for both rules are very similar, both are depicted in Table 1. Note that CANC trivially holds for odd numbers of n, and accordingly verification within the depicted range is achieved in less than 30 seconds.

For the axiom REINF, the run-times can be seen in Nonetheless, the time for finding counter-examples of up to ten voters and five alternatives is significantly lower than for verifying that none exists (see, e.g., the run-times for four and five alternatives in Figure 4, comparing run-times for four and five voters). The smallest counter-examples can be found for both voting rules within less than one minute. These simultaneously serve as proofs, which are both short and easy for a human to inspect.

The smallest example (in number of alternatives) proving that the Black voting rule fails REINF has been found in 48 seconds for three alternatives and five voters:

R 1 = a c c a b b , R 2 = c b a b a c a c b .
The elected alternatives for the profiles R 1 and R 2 are {a, c} and {a, b, c} respectively. For the joined profile R 1 ∪R 2 , the Black voting rule elects the Condorcet winner a instead of the set {a, c} mandated by REINF.

The smallest example, both in number of alternatives and in number of voters, proving that the Copeland rule fails REINF has been found in 32 seconds for three alternatives and four voters:

R 1 = b a a c c b , R 2 = a b b a c c
.

The elected alternatives for the profiles R 1 and R 2 are a and {a, b} respectively. For the joined profile R 1 ∪ R 2 , the Copeland rule elects the set {a, b} instead of only alternative a, which would be required by the REINF property. Both examples can be easily inspected manually.

Literature review

The automatic method presented within this article builds upon the techniques and concepts by [START_REF] Beckert | Automated verification for functional and relational properties of voting rules[END_REF].

It may be used to extend the argumentation framework by [START_REF] Cailloux | Arguing about voting rules[END_REF]. Similarly to our approach, [START_REF] Brandt | Optimal bounds for the no-show paradox via SAT solving[END_REF] apply automated reasoning based on SAT-solvers. They prove existing social choice theorems to obtain minimal bounds for the theorems to hold. More work on efficient applications of automated reasoning in social choice theory exists by [START_REF] Geist | Automated search for impossibility theorems in social choice theory: Ranking sets of objects[END_REF], who automatically detect many impossibility theorems. In another logic-based automatic approach, [START_REF] Xia | Designing social choice mechanisms using machine learning[END_REF] incorporates axiomatic properties to a machine learning framework.

Besides fully automatic approaches, other approaches apply interactive techniques using tactical theorem provers in order to yield unbounded results. This branch was initiated by [START_REF] Nipkow | Social choice theory in hol[END_REF], who provides a mechanized proof for the famous Arrow's possibility theorem. [START_REF] Beckert | Automated verification for functional and relational properties of voting rules[END_REF] use optimized deductive program verification techniques to verify a number of simple voting rules with respect to classical axioms. Moreover, [START_REF] Dawson | Machine-checked reasoning about complex voting schemes using higherorder logic[END_REF] apply interactive theorem proving to formalize and verify complex voting rules according to legal text. Another direction is pursued by [START_REF] Pattinson | Vote counting as mathematical proof[END_REF], who encode complex voting rules in rules of higher-order logic, which serve as proof-carrying code for voting rules and may thus produce certified election results. Verification using tactical theorem provers may lead to very high confidence levels, but requires time-consuming, huge and difficult interactive proofs.

Conclusion and future work

In this article, we have presented and formally defined an automatic approach to argue for and against voting rules based on axiomatic social choice properties. The approach is fully automatic and based on the software analysis technique bounded model checking. Our case study on the voting rules Borda, Copeland and Black shows that illustrative proofs are obtained automatically in reasonable time for small numbers of voters and alternatives. This might appear surprising given the combinatorial structure of preference profiles and set-valued outcomes.

Based on the small-scope hypothesis, we argue that our approach can be used for arguing about voting rules, especially because it provides short and human-readable proofs when it detects that a voting rule fails to satisfy some axiom.

Ideas for future work include extensions to reasoning about rules that are incompletely specified. For example, we could automatically conceive an argument (in the form of an example profile and a set of winners as illustrated in the experiments) that would attack any rule that does not satisfy a given axiom. In this way, it would be possible to also argue against classes of rules, additionally to concrete rules as illustrated in this work. The approach could also be extended to automatically illustrate differences between two sets of axiomatic properties without the need for any explicit voting rule. A more ambitious extension would consist in automatically deriving proofs that a winner is a "right" winner by using several example profiles, as illustrated (using a non-automatic procedure) by [START_REF] Cailloux | Arguing about voting rules[END_REF].

A longer-term objective is to mix the ideas proposed in this work and further optimization techniques [START_REF] Beckert | Automated verification for functional and relational properties of voting rules[END_REF] with elicitation procedures in order to obtain a system that would permit to recommend a voting rule, possibly based on work in elicitation of scoring rules [START_REF] Cailloux | Eliciting a Suitable Voting Rule via Examples[END_REF] and on work that defines justified recommendations as those that resist counter-arguments [START_REF] Cailloux | A formal framework for deliberated judgment[END_REF]].

1

  void dominance(uint prof[N][M], uint win[M]) { 2 uint bad = nondet_uint(), good = nondet_uint(); 3 __CPROVER_assume (0 ≤ bad < M); 4 __CPROVER_assume (0 ≤ good < M); 5 __CPROVER_assume (bad ̸ = good); 6 uint prefergtob[N]; 7 for (uint i = 0; i < N; i++) { 8 uint rankg = M, rankb = M; 9 for (uint r = 0; r < M; r++) { 10 if (prof[i][r] == good) rankg = r; 11 if (prof[i][r] == bad) rankb = r; 12 } 13 prefergtob[i] = (rankg < rankb) ? 1 : 0; 14 } 15 for (uint i = 0; i < N; i++) 16 __CPROVER_assume (prefergtob[i]); 17 __CPROVER_assume (win[bad]); 18 } Listing 1: Pareto-dominance specification for CBMC.

1

  int main(int argc, char * argv[]) { 2 uint prof[N][M], uint win[M];

  denote the "top alternatives" of an elementary profile, meaning, the set of alternatives corresponding to A top . (This is legal as it is unique.) Example 1. The profile R shown below, composed of the linear orders (a, b, c, d) and (b, a, d, c), is an elementary profile corresponding to A top = { a, b }, with A equal to { a, b, c, d }. order ≻ ∈ L (A ), with A of size m, define the cycle corresponding to ≻ as the set of m pairs consisting of the pairs of alternatives (a, b) such that b immediately succeeds to a in ≻ (the rank of b is one more than the rank of a), union the pair of alternatives (z, a) where z and a are respectively the minimal and maximal elements of ≻. For example, the cycle corresponding to the linear order (a, b, c) is {(a, b), (b, c), (c, a)}. Observe that, given any ≻ ∈ L (A ), the cycle corresponding to ≻ corresponds to exactly m linear orders in L (A ). For example, the cycle {(a, b), (b, c), (c, a)} also correspond to (b, c, a) and (c, a, b).

Example 2 .

 2 The profile R shown below is a cyclic profile corresponding to the cycle {(a, b), (b, c), (c, d), (d, a)} with A = { a, b, c, d }.

  for 2, 3, 4 and 5 alternatives in seconds.

Figure 1 :

 1 Figure 1: Run-times for the automatic verification of Pareto-dominance for the Borda voting rule.

  for 2, 3, 4 and 5 alternatives in seconds.

Figure 2 :

 2 Figure 2: Run-times for finding an illustration of the dominance property for the Borda rule.

.

  For the profile R, the Borda rule elects the alternatives {a, c} instead of the Condorcet winner c. Whereas an isomorphic counter-example is found in less than one second for the failure of MAJ, we find two smallest examples (one regarding the number of alternatives, and another one regarding the number of voters) proving the failure of WMAJ. For a minimal number of alternatives, we find the smallest proof for three alternatives and five voters in less than one second:

  for 2, 3, 4 and 5 alternatives in seconds.

Figure 3 :.

 3 Figure 3: Run-times for the verification of REINF for the Copeland voting rule.

  for 2, 3, 4 and 5 alternatives in seconds.

Figure 4 :

 4 Figure 4: Run-times for the verification of REINF for the Black voting rule.

  Figures 3 and 4 for the Copeland rule and the Black rule, respectively. Round dots indicate that violating voting situations have been found for the respective numbers of voters and alternatives. Missing parts (for the case of 5 alternatives) indicate a timeout as defined in Section 5.1.

  An alternative is a Condorcet winner iff it obtains a strict majority against all other alternatives. Any Condorcet winner is unique. COND This property mandates that if there is a Condorcet winner, it becomes the sole winner. MAJ This property requires that, whenever some alternative is placed first by more than half of the voters in (≻ i ) i∈N , it becomes the sole winner. The Condorcet property is stronger than the majority property, i.e., S COND ⊆ S MAJ . Based on MAJ we define a weaker property, the weak majority property WMAJ, for which S MAJ ⊆ S WMAJ .

Table 1 :

 1 Run-times for the verification that Black and Copeland rule each satisfy CANC.

	m	n	verification time
	≤ 3 ≤ 8 4 ≤ 5 and 7 4 6 and 8 5 ≤ 3 and 5 and 7 < 1 minute < 3 minutes < 2 minutes > 30 minutes 5 4 and 6 and 8 > 30 minutes


David González

Francisco Chinesta
email: francisco.chinesta@ensam.eu

Elías Cueto
email: ecueto@unizar.es

D González Aragon

Thermodynamically consistent data-driven computational mechanics

Keywords: Data-driven computational mechanics, GENERIC, Governing Equations

In the paradigm of data-intensive science, automated, unsupervised discovering of governing equations for a given physical phenomenon has attracted a lot of attention in several branches of applied sciences. In this work we propose a method able to avoid the identification of the constitutive equations of complex systems, and rather work in a purely numerical manner by employing experimental data. In sharp contrast to most existing techniques, this method does not rely on the assumption on any particular form for the model (other than some fundamental restrictions placed by classical physics such as the second law of thermodynamics, for instance) nor forces the algorithm to find among a pre-defined set of operators those whose predictions fit best to the available data. Instead, the method is able to identify both the Hamiltonian (conservative) and dissipative parts of the dynamics while satisfying fundamental laws such as energy conservation or positive production of entropy, for instance.

The proposed method is tested against some examples of discrete as well as continuum mechanics, whose accurate results demonstrate the validity of the proposed approach.

Introduction

Scientific experiments in large infrastructures are able to produce petabytes of raw data per day. With the irruption of the Internet of Things and the industry 4.0 paradigm, modern factories are also producing breathtaking amounts of data which should be stored, curated and processed so as to identify tendencies and analyze hidden correlations, patterns and, eventually, make decisions in real time. Therefore, being able to predict equations governing physical phenomena from experimental data has utmost importance. A growing interest has arose in the last years in this direction [START_REF] Bongard | Automated reverse engineering of nonlinear dynamical systems[END_REF] [2] [START_REF] Battista | Numerical investigation of a particle system compared with first and second gradient continua: Deformation and fracture phenomena*[END_REF] [4] [START_REF] Della Corte | Modeling Deformable Bodies Using Discrete Systems with Centroid-Based Propagating Interaction: Fracture and Crack Evolution[END_REF].

There are, however, physical laws of different epistemic nature. One can not attribute the same importance to the second law of thermodynamics than to the response model of a new material, for instance. Often, these last equations are referred to as constitutive equations and have a phenomenological character. Consequently, the process of distilling physical laws from data in an automated way involves the task of finding relevant invariant structures among them -thus satisfying fundamental conservation equations-but also to be able to identify constitutive equations for the different components of the system under scrutiny.

Under non-equilibrium situations the problem becomes even harder. An algorithm able to identify the dynamic structure of the system should be able to identify not only the conserved magnitudes (i.e., to find out its Hamiltonian structure), but also its dissipative nature. Many recent attempts in this sense begin by assuming that the (possibly nonlinear) dynamics of the system is of the form żt = d dt z(t) = f (z(t)), [START_REF] Bongard | Automated reverse engineering of nonlinear dynamical systems[END_REF] where z refers to the pertinent set of state variables, needed for an adequate description of the system. Most existing techniques are designed to search among a set of predefined operators (polynomial terms, sinus, cosinus, ...) those that compone the operator that best fits the experimental results [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] [START_REF] Peherstorfer | Data-driven operator inference for nonintrusive projection-based model reduction[END_REF].

Other techniques, such as [START_REF] Daniels | Automated adaptive inference of phenomenological dynamical models[END_REF], for instance, progressively add terms, somehow resembling a Taylor series, until an adequate degree of accuracy is reached. More recently, the possibility of operating without any constitutive equation has also been explored [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF] [START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF]. In this last work, a set of experimental results substitute constitutive equations by developing a method that operates by finding the experimental result closest to guaranteeing equilibrium. None of the aforementioned techniques, however, impose a priori any fundamental law whose satisfaction is mandatory. Thus, fundamental conservation laws will be satisfied depending on the degree of accuracy of the resulting approximation found by the different algorithms, the noise level in the data, etc. Furthermore, the set of relevant variables z -whose precise meaning will be clear hereafter-could be a priori unknown, or maybe non-relevant variables could have been stored during the experimental campaign. Detecting redundant variables is also an important issue.

In this paper a method is proposed that guarantees the satisfaction of basic principles. Notably, conservation of energy and positive production of entropy are enforced by construction. This is possible by resorting to the so-called GENERIC structure of the problem [START_REF] Oettinger | Beyond Equilibrium Thermodynamics[END_REF]. The GENERIC paradigm generalizes the Hamiltonian structure of classical mechanics by including the dissipative properties of the problem and seems therefore to be the adequate starting point for this project. In Section 2 we review the basics of this paradigm. Then, in Section 3 we introduce a strategy for the identification of the GENERIC constituent blocks from experimental data. In Section 4 we propose three different examples of the performance of the method that reveals the accuracy of the proposed methodology. Finally, the paper is completed with the discussion of the results.

The GENERIC formalism

Under the framework of non-equilibrium thermodynamics, the GENERIC ("General Equation for Non-Equilibrium Reversible-Irreversible Coupling") formalism establishes a completely general equation of the dynamics of a system under reversible and irreversible conditions, as dictated by the evolution of energy and entropy, respectively [START_REF] Grmela | Dynamics and thermodynamics of complex fluids. i. development of a general formalism[END_REF] [START_REF] Oettinger | Beyond Equilibrium Thermodynamics[END_REF]. GENERIC allows for a proper description of the thermodynamical structure of problems within solid and fluid mechanics, thermokinetic effects, chemical reactions or relativistic physics, among other fields in which it has been successfully applied [START_REF] Hans | Nonequilibrium thermodynamics: a powerful tool for scientists and engineers[END_REF]. Models constructed upon the GENERIC formalism are guaranteed to posses a sound thermodynamic structure (they preserve the symmetries of the system and therefore guarantee the subsequent conservation/dissipation laws of energy and entropy, respectively, by Noether's theorem), and therefore it seems to be an appealing choice to start with. In addition, GENERIC paves the way to the development of numerical integration algorithms that can be seen as discrete GENERIC instances, and therefore share these interesting properties [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF] [START_REF] Romero | Algorithms for coupled problems that preserve symmetries and the laws of thermodynamics: Part i: Monolithic integrators and their application to finite strain thermoelasticity[END_REF].

In general, the type of systems we are considering in this work are those that can be described by the laws of classical mechanics. Under this framework, let us consider a closed, isolated system whose evolution is to be described during the time interval I = (0, T ]. The state space of the system will be described by judiciously chosen variables in a space S. The term "judiciously chosen" will be clear hereafter. Therefore, the system will be described at time t ∈ I by a function z t = z(t) :

I → S, z ∈ C 1 (0, T ].
The GENERIC structure of the evolution equations for such a system takes the form żt

= L(z t )∇E(z t ) + M (z t )∇S(z t ), z(0) = z 0 . ( 2
)
L is the so-called Poisson matrix, and will be responsible for the reversible (Hamiltonian) part of the evolution of the system. E represents the energy of the system, as a function of its particular state at time t, z t . In turn, M represents the friction matrix, responsible for the irreversible part of the evolution of the system. Finally, S represents the entropy of the system for the particular choice of variables z. It is important to note that the actual choice of z variables is irrelevant in the sense that any other complete set of variables (possibly including a different number of them) z = f (z) will produce a valid GENERIC structure of the problem, although probably much more complex. Eq. ( 2) is supplemented with the complementary degeneracy conditions, i.e.,

L(z) • ∇S(z) = 0, (3a)
M (z) • ∇E(z) = 0. (3b)
By choosing L skew-symmetric and M symmetric, positive semi-definite, one ensures that

Ė(z) = ∇E(z) • ż = ∇E(z) • L(z)∇E(z) + ∇E(z) • M (z)∇S(z) = 0, (4)
i.e., one ensures the conservation of energy in closed systems.

Equivalently,

Ṡ(z) = ∇s(z) • ż = ∇S(z) • L(z)∇E(z) + ∇S(z) • M (z)∇S(z) ≥ 0, (5)
which ensures the fulfillment of the second principle of thermodynamics.

The reason for choosing GENERIC as the starting paradigm for this work is therefore obvious: it guarantees two basic principles such as the conservation of energy and the second principle of thermodynamics. It is under the umbrella of GENERIC that we will guarantee that the found laws will respect basic principles of utmost importance (such as the symmetries of the system, see [START_REF] Oettinger | Beyond Equilibrium Thermodynamics[END_REF]) and that, after this identification, the resulting time integration schemes will show the right conservation properties [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF].

However, it is important to note that the choice of the state space of the system, z, is not obvious. First, the chosen variables must be measurable experimentally, of course. Otherwise, the identification of the system will be complicated. Secondly, and perhaps more importantly, not every possible description of a system can be cast into the GENERIC framework. For systems obeying the classical laws of mechanics, the most detailed (yet often very impractical) level of description is that of molecular dynamics. By successively coarse-graining the system, keeping in the description those variables that evolve in a slow manifold, one arrives, in the limit, at a description in which the sole variables that are involved are those that are conserved in the system along time, i.e., one arrives at the thermodynamics description of the system. During this journey, different descriptions are possible, and not all are suitable for a GENERIC structure [START_REF] Español | Statistical Mechanics of Coarse-Graining[END_REF]. In general, however, GENERIC seems to be the right framework to express a set of equations yet to be found by means of fitting of experimental results. It involves the most general variables of the system, energy and entropy, while ensuring the satisfaction of the most general requirements of thermodynamics.

However, the proposed methodology, that will be clear in the following section, does not look for the identification of a constitutive equation. Instead, what we propose is to work in a purely numerical, discrete way. In [START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF] we proposed the concept of constitutive manifold, by assuming that the experimental results (typically twelve-dimensional stress-strain couples for solid mechanics or stress-strain rates for fluid mechanics) live in a certain manifold M whose structure is to be identified, possibly with manifold learning techniques. Once identified, the structure of M will allow for an algorithmic numerical approximation of tangent moduli with which to operate in a purely discrete manner, see also [START_REF] Ibañez | Data-driven non-linear elasticity. constitutive manifold construction and problem discretization[END_REF].

The proposed identification algorithm

It is assumed that the available data will provide the necessary information so as to reconstruct the state of the system at discrete time steps. In other words, we assume to know in advance a set of n meas vector-valued measurements of the already mentioned variables of interest Z = {z 0 , z 1 , . . . , z n meas }. By discretizing Eq. ( 2), one arrives at

z n+1 -z n Δt = L(z n+1 )DE(z n+1 ) + M(z n+1 )DS(z n+1 ), (6)
where, for simplicity, we denoted z n+1 = z t+Δt and where L and M are the discrete version of the Poisson and friction operators, respectively. In turn, DE and DS represent the discrete gradients. In general, matrix L is constant over the process, while matrix M frequently varies.

The proposed governing equation identification algorithm will therefore consist of the identification of the discrete values L and M, on one hand, and DE and DS, on the other. In sharp contrast to previous works, such as [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][7], for instance, we do not pursue the identification of a particular term (such as polynomials, sinus, cosinus, ...) in the law. Instead, we will approximate the mentioned ingredients L, M, DE, and DS numerically, and construct with them the already introduced constitutive manifold, see [START_REF] Lopez | A manifold learning approach for integrated computational materials engineering[END_REF][10] [START_REF] Ibañez | Data-driven non-linear elasticity. constitutive manifold construction and problem discretization[END_REF].

Within the realm of variational data assimilation, we will assume that the experimental measurements in Z are exact [START_REF] Manzoni | A reducedorder strategy for solving inverse bayesian shape identification problems in physiological flows[END_REF] [START_REF] Sullivan | Introduction to Uncertainty Quantification[END_REF]. Of course, data will always be corrupted by noise to some extent, and this will oblige us to consider z in a Bayesian framework. This can be done by stating Eq. ( 2) in the framework of a Fokker-Planck equation, which is always possible [START_REF] Soize | The Fokker-Planck equation for stochastic dynamical systems and its explicit steady state solutions[END_REF]. This is currently the focus of our investigations and will be published elsewhere.

It is worth noting that, by identifying the particular form of E(z) we are doing no other thing than identifying what we usually call the "constitutive law" in solid mechanics. In the hyperelasticity framework we usually denote by ψ or W the strain energy density function, and obtain some stress measure (usually the second Piola-Kirchhoff tensor) by means of the computation of the gradient of the elastic energy with respect to some strain measure. In a hyperelastic context we deal with purely conservative systems, and this will be therefore equivalent to obtaining the Poisson matrix and the ∇E(z) terms.

Therefore, the proposed algorithm will consist in solving the following (possibly constrained) minimization problem within a time interval J ⊆ I:

μ * = {L, M, DE, DS} = arg min μ ||z(μ) -z meas ||, (7)
with z meas ⊆ Z, a subset of the total available experimental results. Since L, M, DE, DS will in general be dependent on z, the choice of the size of z meas deserves also some analysis: global fits will in general lead to excessive smoothing of the identified results, while using too small experimental sets will lead to ill-posed minimization problems.

The general form of matrix L is very easy to find, in general -see the examples that follow-while the structure of M is frequently much more intricate. This is why we prefer, in general, to assume L known (if its form can be determined beforehand) and simplify the problem a bit and to re-formulate it as

μ * = {M, DE, DS} = arg min μ ||z(μ) -z meas ||, ( 8
)
so that only the structure of M and the discrete gradients is sought. This avoids to take into consideration all the constraints related to the skew-symmetry of L and improves the convergence of the proposed algorithm. By assuming z approximated in a finite element sense, gradient operators can straightforwardly be cast in matrix form as

DE = Az, (9a) DS = Bz, ( 9b
)
where A and B represent the discrete, matrix form of the gradient operators leading to DE and DS, respectively. Once identified, the manifold structure of {L, M, DE, DS} will also serve to construct, following Eq. ( 6), energy-conserving integration schemes with non-decreasing entropy, see [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF].

At this point, several identification strategies can be envisaged. These are summarized next.

Identification strategies

To solve problem [START_REF] Daniels | Automated adaptive inference of phenomenological dynamical models[END_REF], several strategies could be implemented, depending on the usage of monolithic or staggered approaches, constrained or not.

Monolithic vs. staggered approaches

For some systems, it is possible to easily design experiments under conservative (purely Hamiltonian) conditions. For instance, some systems can be isolated so as to obtain isothermal results with whom to identify the Poisson term of the GENERIC form of the model. Subsequent non-equilibrium experiments will help to identify the dissipative structure of the system. Of course, this is not always possible. One of such examples is that of viscoelastic materials. Viscous and elastic (purely conservative) parts of their behavior can not be decoupled experimentally. In those cases, the whole set of ingredients of the GENERIC formalism must be identified at once.

Total vs. partial identification

We have already mentioned that, although Eq.( 8) incorporates the four ingredients of GENERIC to be identified (i.e., L, M, DE, DS), some of these matrices are easily identified without resorting to any experimental result. For instance, if the system is of mechanical nature and purely Hamiltonian, by choosing z = {q(t), p(t)} (position and momentum of each degree of freedom), the Poisson matrix takes always the form

L(z) = 0 1 -1 0 . ( 10
)
This is the case very often. In other situations, very simple considerations shed light on the form of L and -less frequently-of M .

Constrained or unconstrained minimization

As mentioned before, one of the salient features of GENERIC is the guaranteed conservation of energy and positive entropy production, see Eqs. ( 4)- [START_REF] Della Corte | Modeling Deformable Bodies Using Discrete Systems with Centroid-Based Propagating Interaction: Fracture and Crack Evolution[END_REF].

In order to identify the form of the discrete gradient operators, DE, DS, two possibilities arise: to verify a posteriori that the numerically identified terms fulfill these requirements, or to impose it a priori, by solving a constrained minimization procedure:

μ * = {M, DE, DS} = arg min μ ||z(μ) -z meas ||, (11a) subject to: L • Bz = 0, (11b)
M • Az = 0, (11c)
where restrictions (11b) and (11c) will help, in general, to facilitate the identification procedure while ensuring the thermodynamic consistency of the resulting model.

Results

To validate the proposed methodology, we begin by considering two different finite-dimensional problems whose GENERIC structure is known, so as to be able to determine its degree of accuracy. This allows us to determine the best strategies in order to minimize errors in the identification procedure. Then, an infinite-dimensional problem is analyzed that also proves the validity of the suggested strategy.

Simple pendulum

The first problem corresponds to a simple pendulum whose purely Hamiltonian structure, in absence of any dissipation, responds to Eq. [START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF]. Therefore, only the discrete matrix form of the gradient, A, was identified. The accuracy of the proposed identification method was tested by reproducing the dynamics of the system, integrated by constructing the time integrators proposed in [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF], namely, the simple scheme outlined by Eq. ( 6). As can be noticed, the agreement is almost perfect. Knowing in advance the form of L, the identification of A employing all the available data (241 time increments) lead to an error vanishing to machine precision, see Fig. 1.

On the contrary, if the identification procedure is applied by taking subsets of five experimental measurements, an error in L2-norm of 0.049% was obtained, see Fig. 2.

Thermoelastic pendulum

As an example of a non-equilibrium problem to start with, we have chosen a problem presented in [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF]. The system consists of two pendula connected by thermoelastic springs. Two masses m 1 and m 2 are thus connected by springs of internal energy e 1 and e 2 , and are allowed to oscillate around a fixed point, see Fig. 3. We employ the classical notation of q i , p i , i = 1, 2 for position and momenta, respectively. In turn, springs have entropies s j , and longitudes at rest λ 0 j , j = a, b. Under this perspective, the state space of the system can be defined as Fig. 2 Results for the pendulum problem. Data (circles) vs. numerical integration of the identified system when the identification procedure is applied every five time increments. The GENERIC structure of the problem will naturally involve the internal energy of the system, as already discussed. It can be computed as the sum of the kinetic energy (coming from the movement of the masses) and the potential energy in the springs, i.e.,

S = {z = (q 1 , q 2 , p 1 , p 2 , s 1 , s 2 ) ∈ (R 2 ×R 2 ×R 2 ×R 2 ×R×R), q 1 = 0, q 2 = q 1 }.
E(z) = K 1 (z) + K 2 (z) + e a (λ a , s a ) + e b (λ b , s b ),
where

λ a = √ q 1 • q 1 , λ b = (q 2 -q 1 ) • (q 2 -q 1 ).
It is assumed that the temperature in the springs, θ j , changes due to the Joule effect, giving rise to

θ j = ∂e j ∂s j , j = a, b.
Assuming a conductivity κ in the springs, the resulting equations of the system will be qi =

p i m i , ṗi = - ∂ ∂q i (e a + e b ), ṡj = κ θ k θ j -1 , with i = 1, 2, j = a, b, k = j.
This gives rise to the following gradients, that are to be identified,

∇E(z) = f a n a -f b n b , f b n b , p 1 m 1 , p 2 m 2 , θ a , θ b ,
∇S(z) = (0, 0, 0, 0, 1, 1), with f j , n j , j = a, b, the forces in the springs and their respective unit vector along their direction. Note that we assume them unknown and employ these expressions as ground truth just to test the accuracy of the proposed methodology.

Something which is particularly appealing from the structure of GENERIC, already commented in the introduction, is that different choices of state-space variables are equally valid. For instance, one can imagine S = {z | z = (q 1 , q 2 , p 1 , p 2 , e 1 , e 2 )}, or even S = {z | z = (q 1 , q 2 , p 1 , p 2 , θ 1 , θ 2 )}, and it will produce equivalent results (although with different complexity in the resulting expressions, see [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF]). Under the first of the choices, Poisson and friction matrices will result in

L(z) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
0 0 1 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , M (z) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

κ θ b θ a -κ 0 0 0 0 -κ κ θ a θ b ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
, although, as explained before, we assume M unknown.

The different identification procedures suggested in Section 3 allowed to determine the precise form of the gradients of energy and entropy, respectively. Errors in the numerical reconstruction of the pendula trajectories with the just identified terms are reported in Table 1. In that table, the mentioned number of points refers to the size of data sets employed for the identification, i.e., the size of the set z meas in Eq. ( 8).

An example of the reconstructed trajectory of the masses is shown in Fig. 4. Fig. 4 Results for the thermal pendulum problem. Trajectories of the two masses. Data (circles) vs. numerical integration of the identified system. Monolithic identification on sets of |z meas | = 34 data points.

Oldroyd-B fluid

The field of non-Newtonian fluid mechanics is where GENERIC originated and one of the most active in terms of identification of complex constitutive equations. The Oldroyd-B fluid model is a clear example of an equation that can be derived by starting from very different assumptions: from continuum considerations to coarse-graining microscopical descriptions or, finally, by employing a Fokker-Planck probabilistic description [START_REF] Owens | Computational Rheology[END_REF].

The Oldroyd-B model originated as a means to model the flow of diluted polymeric solutions. These polymer chains can be considered both from a purely macroscopic point of view as well as from a microscopic one, by modeling them as linear dumbbells in a Newtonian medium. From a purely continuum point of view, the constitutive equation of an Oldroyd-B fluid arises from considering the deviatoric part T of the stress tensor σ (the so-called extra-stress tensor), of the from

T + λ 1 ∇ T = η 0 γ + λ 2 ∇ γ , ( 12
)
where the triangle denotes the non-linear upper-convected derivative introduced by Oldroyd [START_REF] Walters | The distinctive CFD challenges of computational rheology[END_REF]. The coefficients η 0 , λ 1 and λ 2 are material constants.

Finally, following standard notation, γ = (∇ s v) = D represents the strain rate tensor.

If we now consider separately the stress in the solvent (denoted by a subscript s) and polymer (denoted by a subscript p) components as

T = η s γ + τ , this gives rise to τ + λ 1 ∇ τ = η p γ,
which is the constitutive equation for the elastic stress.

An alternative derivation can be obtained by starting from a population of linear dumbbells immersed in a Newtonian fluid [START_REF] Owens | Computational Rheology[END_REF]. What is here intended is to take the Oldroyd-B model as a paradigm of a class of viscoelastic (and therefore non Hamiltonian) fluids that encompasses both macroscopic and mesoscopic models. Therefore, we assume that we have a set of pseudoexperimental measurements arising from an Oldroyd-B fluid and try first to identify its constitutive equation in a GENERIC framework, and second, to reconstruct the results by integrating the just found constitutive equations by means of Eq. ( 6).

GENERIC form of the problem

In a GENERIC framework, rather than dealing with constitutive equations for the (elastic) extra-stress tensor, as one intuitively usually starts by, it is often suggested (see, for instance, [START_REF] Pasquali | Theoretical modeling of microstructured liquids: a simple thermodynamic approach[END_REF] and references therein) to start by using the conformation tensor c = rr , as the second moment of the dumbbell endto-end distance distribution function. In fact, this choice allows for a natural introduction of thermal fluctuations in the model [START_REF] Vázquez-Quesada | Consistent scaling of thermal fluctuations in smoothed dissipative particle dynamics[END_REF]. It can be demonstrated that for modeling viscoelastic fluids, the only variables needed to describe the state of the system are the momentum density p = ρv, and the conformation density, ĉ = ρc p c, where v represents the fluid's velocity, ρ is the density and c p is the number of polymer molecules per unit mass, see [START_REF] Vlasis | Atomistic monte carlo simulations of polymer melt elasticity: their nonequilibrium thermodynamics generic formulation in a generalized canonical ensemble[END_REF] and references therein. This choice can result somewhat unnatural for scientists used to continuum models. And indeed it is, if one takes into account that the conformation tensor is not easily observable from an experimental point of view. However, the reader can easily notice its similitude to the left Cauchy-Green or Finger strain tensor b = F • F T . Indeed, in a crosslinked rubber, c = N 2 b/3, where N is the number of statistical segments in a polymer chain and its length [START_REF] Pasquali | Theoretical modeling of microstructured liquids: a simple thermodynamic approach[END_REF].

It is well known that the extra-stress tensor takes the form

τ = - n 2ζ 12 ∇ c,
i.e., it is in fact proportional to the convective (or Oldroyd) derivative of the conformation tensor. This magnitude is more suitable as a characteristic internal variable of the fluid than c = rr itself.

Therefore, in order to model a viscoelastic fluid in the GENERIC framework, we must consider, for each of the M fluid particles that compose the (discretized) fluid, their position r j , velocity v j , energy E j and extra-stress tensor τ j . The internal energy E j includes, of course, the kinetic energy contribution from the solvent and the beads, plus the potential energy of solvent-solvent, solvent-bead and bead-bead interaction. Each one of the M fluid particles is assumed to contain j = 1, 2, . . . , N dumbbells and, of course, an entropy S j (E j , τ j , V j ), with V j the volume associated to each particle.

Therefore, using the previous notation, the set of variables employed to describe the phase space of the fluid will be

S = {z = (r j , v j , E j , τ j , j = 1, 2, . . . , M) ∈ (R 2 × R 2 × R × R 4 ) M }. (13)

Implementation details

Although the just introduced strategy is entirely equivalent to the existing data-driven simulation approaches by Ortiz and coworkers [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF] [28], in the sense that it performs a constrained minimization problem at every integration point, alternative approaches exist. Instead of fitting point-wise the best possible constitutive values provided by data, it is possible to do it in a non-local manner, by employing model order reduction approaches.

Model reduction consists, essentially, in finding the slow manifold in which the system evolves. Therefore, for a (full order) problem formulated under the GENERIC paradigm, ensuring that its reduced-order version also verifies the GENERIC structure is not obvious. The interested reader can consult [START_REF] Christian | Preservation of thermodynamic structure in model reduction[END_REF]. The obvious advantage is the drastic reduction in the number of degrees of freedom.

In this case we employed a very simple strategy based on the employ of the Proper Orthogonal Decomposition (POD) of the pseudo-experimental results [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF] [31] [START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF]. An alternative, perhaps more accurate, is to employ some form of non-linear dimensionality reduction, or manifold learning strategy, as in [START_REF] Millán | Nonlinear manifold learning for model reduction in finite elastodynamics[END_REF].

Our pseudo-experimental results were generated with a CONNFFESSIT approach [START_REF] Laso | Calculation of viscoelastic flow using molecular models: the connffessit approach[END_REF] [START_REF] Elias Cueto | Meshless stochastic simulation of micro macro kinetic theory models[END_REF] to the problem of a Couette flow of an Oldroyd-B fluid. The Couette flow is the pure shear flow between two rigid plates produced typically in rheometers, and thus the most frequent form of experimental results for this type of fluids. The CONNFFESSIT approach combines finite elements and stochastic techniques for the solution of micro-macro flows of viscoelastic flows.

Given the rheometric character of the flow (in other words, its homogeneous character), the code discretizes only the one-dimensional distance h between the rheometer plates, by employing one hundred finite elements (and, consequently, 101 nodes). We consider the startup of a flow with a Reynolds number of 0.1, and Weissemberg number 0.5. Each finite element node is equipped with 10000 dumbbells.

A simple Proper Orthogonal Decomposition (POD) of the 150 snapshots of the velocity profile, see Fig. 5, revealed that the flow can be expressed in terms of a few degrees of freedom, see Fig. 6. In fact, only six modes produce a decay on the eigenvalues of more than four orders of magnitude. Therefore, instead of the 101 nodes employed to build the mesh, our identification procedure employed the just mentioned six degrees of freedom. Thus, for each of the six retained degrees of freedom, the Poisson and friction matrices are known to have the form

L k (z) = ⎛ ⎜ ⎜ ⎝ 0 1 0 0 -1 0 1 -1 0 -1 0 0 0 1 0 0 ⎞ ⎟ ⎟ ⎠ , M k (z) = ⎛ ⎜ ⎜ ⎝ 0 0 0 0 0 1 1 0 0 1 1 0 0 0 0 1 ⎞ ⎟ ⎟ ⎠ , k = 1, . . . , n POD ,
with n POD the number of employed POD modes. They are, however, assumed unknown and their form is determined by the proposed method. In this case, the identification includes both L and M without any significant difficulty.

Note that, due to the particular form of the Couette flow, the variables in S, see Eq. ( 13), turn out to be one-dimensional, and hence the 4 × 4 size of the Poisson and friction matrices.

With these ingredients, the pseudo-experimental and identified results for nodes 1, 11, . . . , 101 along the vertical axis of the flow, are represented in Fig. 7.

As in previous examples, the accuracy of the proposed methods depends on the number of time samples employed to perform the minimization problem. In this case, by performing a constrained minimization (thus enforcing conservation of energy and strictly positive entropy dissipation) with sets of 150 snapshots, the error in the velocity field turned out to be of 10 -8 . The unconstrained minimization strategy, on the contrary, was unable to find any suitable solution for the problem.

Numerical integration of an arbitrary flow

Of course, the final goal of this procedure is not to reproduce the experimental results but to be able to integrate different flow regimes once the constitutive manifold of the particular fluid has been identified. Assume, for simplicity, that we have identified the discrete gradient operators at time t n , DE = A n z n and DS = B n z n , and the Poisson and friction matrices L and M, respectivelythat turned out to be constant for this problem-for two flows under different conditions, say Weissemberg, W e, numbers. We have, therefore, the values L, M and the triplets A 1 n , B 1 n , z 1 n and A 2 n , B 2 n , z 2 n , together with the current state of the system, z n . The objective is then to numerically integrate, without any explicit knowledge on the constitutive equations of an Oldroyd-B fluid, the state equations for the flow under a third W e value, characterized by A obj , B obj . Assume known the value of the weights w 1 and w 2 that will make Eq. ( 6) hold under the new flow conditions, namely,

z obj n+1 = z n + Δt • LA obj + M B obj z n = z n + w 1 Δt • LA 1 n + M B 1 n z 1 n + z n + w 2 Δt • LA 2 n + M B 2 n z 2 n
with L and M also known. Some simple algebra leads to the following conditions:

A obj z n = w 1 A 1 n z 1 n + w 2 A 2 n z 2 n , B obj z n = w 1 B 1 n z 1 n + w 2 B 2 n z 2 n .
Computing the interpolating weights w 1 , w 2 , . . . , w , with the number of neighbouring experimental results to interpolate with, can be done, for instance, by resorting to non-linear dimensionality reduction techniques such as Locally Linear Embedding (LLE), which produces these weights as a byproduct of the dimensionality reduction process [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]. By doing so, the results of integrating a Couette flow at W e = 0.75 from the experimental results obtained at W e = 0.5 and W e = 1.0 (thus, w 1 = w 2 = 0.5 in a straightforward manner, for this simple illustration) are shown in Fig. 8. In it, the data-driven results (we recall again that no constitutive equation was employed to perform the integration) are compared towards the result of a direct numerical integration. Results for the stress component τ xy are shown. L2-norm error of the velocity field showed to be 0.4988% in this particular case.

Poiseuille flow

Obviously, the aim of the just developed technique is not to be able to reproduce experimental results. Rather, it is to be able to simulate situations for which no experimental result exist. Usually, rheologists employ rheometers to fully characterize the constitutive equations of a particular fluid. The simplest of these rheometers provokes simple shear (Couette) flows, as in the previous sections. In this section we will try to demonstrate that, once the constitutive manifold has been properly characterize by means of the necessary experimental results, the approach just developed could be employed for the simulation of flows for which no previous result exists. So, let us assume that we have performed enough experiments so as to characterize the constitutive manifold. In this case, we have simulated a big amount of shear (Couette) flows to characterize the fluid. With the help of these (pseudo-)experimental results, we are about to simulate the startup of a Poiseuille flow. By taking a time increment of Δt = 0.006, we interpolate from the set of available (pseudo-)experimental results the A and B matrices that are employed to advance the flow in time. Figure 9 represents the evolution of the τ xy component along time, both for the reference solution and the just computed one. For this flow, the mean relative error for this value is between 4.13% and 16.14% for the different nodes of the discretization (which is actually one-dimensional). Of course, these errors can be lowered by employing finer meshes and smaller time steps.

Discussion

Within the realm of data-intensive computational mechanics, this paper presents a novel methodology able to integrate the equations of motion without resorting to any constitutive equation (the so-called data-driven computational mechanics framework first established in [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF]). The main novelty in the approach herein consist in the development of a thermodynamically consistent approach to the problem that guarantees the fulfillment of the first and second principles of thermodynamics, independently of the quality of the experimental results. This is in sharp contrast to existing approaches in the literature, either the ones working without any constitutive assumption, see [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF] [18] [10] [START_REF] Ibañez | Data-driven non-linear elasticity. constitutive manifold construction and problem discretization[END_REF], or those that identify the form of the governing equations by choosing among a predefined set of equations, [6] [7]. The present approach is based on the usage of the GENERIC paradigm, a framework that guarantees the strict energy conservation and positive entropy production, along with the symmetries of the system, and that generalizes the Hamiltonian form of classical mechanics to non-equilibrium systems. The proposed methodology thus identifies the proper GENERIC structure of the problem. The method identifies either the whole GENERIC structure composed by the Poisson and friction matrices of the problem, plus the precise numerical form of the gradients of the energy and entropy. However, for some complex problems it has been found that the form of the Hamiltonian and dissipative structures of the problem is easy to identify, and that the method works considerably better by just finding the form of the gradients of energy and entropy, respectively. This is equivalent, of course, to finding in a numerical, point-wise manner, the form of the constitutive manifold that the authors defined in a previous work, see [START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF].

But the proposed methodology not only provides with a suitable method to unveil the precise form of the Hamiltonian and dissipative structures of the problem. It also provides a robust and thermodynamically consistent way to integrate these equations without employing any constitutive equation (the so-called data-driven computational mechanics paradigm).

The aforementioned properties have been demonstrated by solving different problems: from simple pendulums with discrete degrees of freedom to complex viscoelastic fluids. In this last case, the proposed technique has been combined with a classical model order reduction technique (POD) that allows to highly simplify the problem and to greatly diminish the number of degrees of freedom and variables to identify.

Fig. 1

 1 Fig.1Results for the pendulum problem. Data (circles) vs. numerical integration of the identified system.

Fig. 3

 3 Fig. 3 Double thermal pendulum.

Fig. 5

 5 Fig.[START_REF] Della Corte | Modeling Deformable Bodies Using Discrete Systems with Centroid-Based Propagating Interaction: Fracture and Crack Evolution[END_REF] Velocity profile at every Δt = 0.01 seconds (t ∈ (0, 1.5]) for the startup of a Couette flow (shear flow between two plates, typical of a rheometer). In this case, velocity is imposed at the bottom plate, with h = 1.0 and v shear = 1.0. Notice the overshoot of the velocity profile just before the permanent regime (linear velocity profile).

6

 6 POD modes and related eigenvalues for the Couette flow.

Fig. 7

 7 Fig. 7 Experimental vs. reconstructed phase space variables for the Oldroyd-B fluid. (a) Displacement, (b) velocity, (c) Energy, and (d) xy component of the extra-stress tensor, τ xy . Every variable is represented for nodes 1, 11, . . . , 101 along the height of the channel, h.

Fig. 8

 8 Fig. 8 Stress tensor component τ xy evolution for the reconstructed flow versus reference solution. Results for the points at y = 0, y = h 2 and y = h are shown.

Fig. 9

 9 Fig. 9 Evolution of the τ xy component for a given point in a Poiseuille flow. Red asterisks, reference values. Blue circles, obtained values.

Table 1

 1 Errors (L2-norm) for the double thermal pendulum.

			Strategy
		Staggered	Monolithic
		84 points	17 points	167 points	34 points
	Constrained	0.000001%	0.0684%	0.0023%	0.3072%
	Unconstrained 0.000018%	0.3114%	0.00036%	0.0399%

Acknowledgements This work has been supported by the Spanish Ministry of Economy and Competitiveness through Grants number DPI2017-85139-C2-1-R and DPI2015-72365-EXP and by the Regional Government of Aragon and the European Social Fund, research group T88.


D González

F Chinesta

E Cueto
email: ecueto@unizar.es

Consistent Data-Driven Computational Mechanics

We present a novel method, within the realm of data-driven computational mechanics, to obtain reliable and thermodynamically sound simulation from experimental data. We thus avoid the need to fit any phenomenological model in the construction of the simulation model. This kind of techniques opens unprecedented possibilities in the framework of data-driven application systems and, particularly, in the paradigm of industry 4.0.

INTRODUCTION

In recent years an increasing interest on the development of methods for data-driven simulation has arose. These methods eliminate the need for phenomenological constitutive equation fitting, a process that is frequently cumbersome. Instead, they work operating solely on experimental data [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF][START_REF] Ibañez | Data-driven non-linear elastic-ity. constitutive manifold construction and problem discretization[END_REF][START_REF] Ibanñez | A manifold learn-ing approach to data-driven computational elasticity and inelasticity[END_REF][START_REF] Lopez | A manifold learning approach for inte-grated computational materials engineering[END_REF]. However, the need to guarantee thermodynamic consistency continues to be an issue for these methods. Since in general any experimental result will incorporate errors or some level of experimental noise, a doubt is frequently expressed in terms of the need to guarantee the satisfaction of basic principles. These include the laws of thermodynamics, notably.

To overcome these limitations, in a recent work, we have proposed to formulate the problem in the so-called GENERIC formalism [START_REF] Grmela | Dynamics and thermo-dynamics of complex fluids. i. development of a general formalism[END_REF][START_REF] Oettinger | Beyond Equilibrium Thermodynamics[END_REF][START_REF] Hans | Non equilibrium thermodynamics: a powerful tool for scientists and engineers[END_REF]. The GENERIC ("General Equation for Non-Equilibrium Reversible-Irreversible Coupling") formalism establishes a completely general equation of the dynamics of a system under reversible and irreversible conditions, as dictated by the evolution of energy and entropy, respectively.

The state space of the system will be described by judiciously chosen variables in a space S. Therefore, the system will be described at time t ∈ I by a function z t = z(t) : I → S, z ∈ C 1 (0,T). The GENERIC structure of the evolution equations for such a system takes the form This equation must be supplemented with the so-called "degeneracy conditions", By choosing L skew-symmetric and M symmetric, positive semi-definite, one also ensures the conservation of energy and the strict entropy dissipation.

PROPOSED IDENTIFICATION ALGORITHM

We assume to know a certain amount of experimental data in the form of a vector-valued set Z={z 0 , z 1 , …, z nmeas }. By discretizing the GENERIC equation, one arrives at where L, and M represent the discrete version of the Poisson and friction operators, L and M, and DE and DS represent the discrete gradients, respectively. It is worth noting that, by identifying the particular form of E(z) we are doing no other thing than identifying what we usually call the "constitutive law" in solid mechanics.

Therefore, in the just proposed data-driven identification algorithm we are solving a data assimilation problem of the form with z meas a subset of Z, the total number of experimental results.

Simple pendulum

The first problem corresponds to a simple pendulum whose purely Hamiltonian structure, in absence of any dissipation, responds to the case M=0. Therefore, only the discrete matrix form of the gradient DE was identified. In this case, the error in the identification procedure vanishes within machine precision, if we employ the whole set of 241 experimental measurements. By taking them in sub-sets of five experimental values, the error increases to 0.44%, see Fig. 2.

Thermoelastic pendulum

As an example of non-equilibrium system, we have chosen the thermoelastic pendulum introduced in [START_REF] Romero | Thermodynamically consistent time-stepping algo-rithms for non-linear thermomechanical systems[END_REF], see Fig. 3. The set of variables describing the pendulum will be given by

S = {z = (q 1 , q 2 , p 1 , p 2 , s 1 , s 2 )∈(ℜ 2 ×ℜ 2 ×ℜ 2 ×ℜ 2 ×ℜ×ℜ), q 1 ≠0, q 2 ≠0},
Where q and p represent position and momenta, respectively, and s represents entropy. Results for the identification of this system are shown in Fig. 4.

CONCLUSIONS

The present approach is based on the usage of the GENERIC paradigm, a framework that guarantees the strict energy conservation and positive entropy production, along with the symmetries of the system, and that generalizes the Hamiltonian form of classical mechanics to non-equilibrium systems. The proposed methodology thus identifies the proper GENERIC structure of the problem. The method identifies either the whole GENERIC structure composed by the Poisson and friction matrices of the problem, plus the precise numerical form of the gradients of the energy and entropy. However, for some complex problems it has been found that the form of the Hamiltonian and dissipative structures of the problem is easy to identify, and that the method works considerably better by just finding the form of the gradients of energy and entropy, respectively. This is equivalent, of course, to finding in a numerical, point-wise manner, the form of the constitutive manifold that the authors defined in a previous work.

But the proposed methodology not only provides with a suitable method to unveil the precise form of the Hamiltonian and dissipative structures of the problem. It also provides a robust and thermodynamically consistent way to integrate these equations without employing any constitutive equation (the so-called data-driven computational mechanics paradigm).

FIGURE 1 .

 1 FIGURE 1. Results for the identification of the Hamiltonian structure of a simple pendulum.

FIGURE 2 .

 2 FIGURE 2. Results for the identification of the Hamiltonian structure of a simple pendulum. Data sets composed on only five temporal series.

FIGURE 3 .

 3 FIGURE 3. Thermoelastic pendulum.

FIGURE 4 .

 4 FIGURE 4. Results for the thermal pendulum problem. Experimental results were grouped in sets of 34 time instants.

ACKNOWLEDGMENTS

This work has been supported by the Spanish Ministry of Economy and Com-petitiveness through Grants number DPI2017-85139-C2-1-R and DPI2015-72365-EXP and by the Regional Government of Aragon and the European Social Fund, research group T88.