# Data-driven material modeling based on the Constitutive Relation Error

Pierre Ladevèze, Ludovic Chamoin

► To cite this version:

Pierre Ladevèze, Ludovic Chamoin. Data-driven material modeling based on the Constitutive Relation Error. Advanced Modeling and Simulation in Engineering Sciences, 2024, 11 (1), pp.23. 10.1186/s40323-024-00279-x. hal-04847388

HAL Id: hal-04847388

https://hal.science/hal-04847388v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Data-driven material modeling based on the Constitutive Relation Error

Pierre Ladevèze<fnref n="1" /> and Ludovic Chamoin<fnref n="2" />,<fnref n="1" />*

**Abstract**

Prior to any numerical development, the paper objective is to answer first to a fundamental question: what is the mathematical form of the most general data-driven constitutive model for stable materials, taking maximum account of knowledge from physics and materials science? Here we restrict ourselves to elasto-(visco-)plastic materials under the small displacement assumption. The experimental data consists of full-field measurements from a family of tested mechanical structures. In this framework, a general data-driven approach is proposed to learn the constitutive model (in terms of thermodynamic potentials) from data. A key element that defines the proposed data-driven approach is a tool: the Constitutive Relation Error (CRE); the data-driven model is then the minimizer of the CRE. A notable aspect of this procedure is that it leads to quasi-explicit formulations of the optimal constitutive model. Eventually, a modified Constitutive Relation Error is introduced to take measurement noise into account.

**Keywords:** Data-driven modeling, Materials science, Constitutive Relation Error, Elasto-(visco-)plasticity

**Introduction**

A first pillar in computational mechanics is a certain abstract representation and understanding of the world we live in, in terms of models described by physics laws but also sometimes debatable theories and knowledge developed over the centuries. In this framework, the equilibrium equations, the compatibility equations and the two principles of the thermodynamics of irreversible processes constitute a foundation on which, in our opinion, any approach should be based. A second pillar is experimental data and its impact is continuously increasing nowadays. Data can be of different natures; indeed, it may refer to the geometry of the studied structure, its loading, or the response of the materials constituting the structure. Any calculation model is situated between these two pillars with more or less weight for one or the other. When the influence of data is largely dominant, and no questionable hypothesis is added in the modeling, we generally speak of data-driven approaches. This way, the use of parameterized behavior models from materials science does not seem to us to be a data-driven approach.

Here, we precisely focus on the material. We consider the data-driven modeling of stable complex materials which has become a growing trend. The pioneer works were due to M.

<figure><img src="image_3.png" /></figure>

<footnote>*Correspondence:
ludovic.chamoin@ens-paris-saclay.fr</footnote>

<footnote n="1">¹CentraleSupélec, ENS
Paris-Saclay, CNRS,
LMPS-Laboratoire de Mécanique
Paris-Saclay, Université
Paris-Saclay, 4 Avenue des
Sciences, 91190 Gif-sur-Yvette,
France</footnote>

<footnote n="2">²Institut Universitaire de France
(IUF), 1 rue Descartes, 75231
Paris, France</footnote>
Ortiz and F. Chinesta and their co-workers [39, 41, 42, 53, 57]. In these works, materials science was either discarded or understated. Further works following this line introduced a structuration of data [15]. Alternative developments have recently emerged, in which some physics background is involved when learning the constitutive model. Pioneer works in this direction are [26, 40] and the research topic is evolving fast with various techniques coupling deep learning and physics information to describe the behavior of materials using neural networks e.g., [8, 28, 38, 44, 55, 68, 76, 78].

Looking deeper in the bibliography of “physics-augmented” learning, in which physics knowledge (e.g., principles of thermodynamics or symmetries) is integrated in the neural network architecture, a large list of works is available e.g., [9, 31, 32, 54] to cite some very recent ones. We may separate contributions depending on the data used. A first set of contributions aims to train neural networks in a supervised learning procedure with a strain-stress (or strain-free energy) database generated from a known constitutive model [3, 23, 60, 76]; nevertheless, getting such measurement couples under complex loading is hardly possible in practice. Another set of contributions addresses unsupervised training of neural networks for constitutive modeling [8, 38, 75]; for instance, in the EUCLID approach proposed in [21, 75], full-field displacement observations are processed to derive strain and stress fields employed as the learning input with a loss function that penalizes the non-satisfaction of equilibrium.

In contrast to the previous pioneer approaches, the aim of this paper is to define a data-driven mathematical model prior to any numerical development, taking maximum account of knowledge on physics and materials science. This paper is thus limited to some fundamental aspects; numerical treatment, which involves solving optimization, regression and interpolation problems for which machine learning tools (and in particular deep neural networks) complement the more conventional tools, is not covered. Here we restrict ourselves to stable elasto-(vico)plastic materials under the small displacement assumption. The central concept of the proposed approach is a tool, referred to as the Constitutive Relation Error (CRE), the computed data-driven constitutive model being the minimizer of the CRE.

The CRE has been initially introduced in model verification for linear problems [45]. Precisely, the CRE is built on the separation of the governing equations of the mechanical or physical problem, defined over the space-time domain: equilibrium and compatibility equations on the one hand, and constitutive equations on the other hand. In contrast to the constitutive relations which are empirical and strongly depend on experimental data, the other equations can be qualified as “exact”. Suitable approximations should thus be constructed that “exactly” satisfy the “exact equations”; their quality can be assessed by measuring the mismatch in constitutive equations by the CRE tool. Among the numerous works on CRE for model verification, one can mention the reviews available in the book [48], and the chapters [12, 49]. Direct applications of the CRE for model updating in elasticity and elastodynamics have been performed in the 1990s, in the context of parametrized constitutive relations; one can cite [6, 10, 22, 25, 43]. An extension of the CRE to unreliable data by using a modified formulation (mCRE) has been introduced in [13, 46] for the updating of vibration models. Since then, numerous further works have been performed over the years for model updating and validation for linear or nonlinear materials, in quasi-statics but also in dynamics e.g., [2, 7, 11, 16, 17, 19, 37, 59, 66, 72] to cite a few of them. Recent applications to sequential data assimilation and structural health
monitoring by complementing mCRE with Kalman filtering have also been investigated lately [18,58].

Like all inverse problems, data-driven material modeling is characterized by an over-abundance of data and therefore leads to the solution of ill-posed problems. CRE was primarily used to provide a solution to such problems (in fact to the usual family of ill-posed problems called standard), in situations where the constitutive relation was known; the goal was then to identify or update the associated parameters. The experimental data derived from full-field measurements (even though sparser measurements can also be considered) is related to the family of structures tested, made from the same material and subjected to different loads. Here, we show that the CRE framework can again be used to calculate the data-driven model that is defined for stable materials by thermodynamic potentials, the calculated potentials being the CRE minimizer. The advantage of the approach is the use of an error metric which is directly focused on what needs to be learnt, that is the constitutive relation. To start, we consider linear and nonlinear elastic materials with convex energy, for which the CRE approach gives original solutions. This approach is compared with that proposed in the pioneer work [41].

For elasto-(visco-)plastic materials, materials science is a great help. The very common internal variables framework may be advantageously used but data-driven approaches involve two main difficulties. First, the hidden internal state variables are not known a priori, nor is their number. Moreover, they are not intrinsic, and transformations on these variables leading to the same material model are possible. This is why, for stable materials with a convex energy, we introduce the so-called normal formulation of the hidden internal variables resulting from a transformation of these variables while satisfying, of course, the two principles of the thermodynamics of irreversible processes. In addition, we follow what is classically done in materials science to describe stable elasto-(visco-)plastic materials by referring to the sub-family of standard materials introduced in [29]. It follows that the mathematical shape of data-driven models for stable elasto-(visco-)plastic materials is defined by only one convex function. This formulation, both normal and standard, is at the core of the proposed data-driven CRE approach. A unique solution is given to the standard ill-posed problems, as the minimizer of the CRE functional; it is defined by its own constitutive relation, which differs from the reference relation. The optimal data-driven model is then computed by a constrained minimization of the CRE using convex optimization.

Eventually, we consider noisy data. Following [33,71], we assume that the main source of noise comes from full-field measurements. A modified Constitutive Relation Error (mCRE) is then introduced, still for elasto-(visco-)plastic materials, and developments are very similar to the case of noiseless data for the derivation of a data-driven constitutive model.

To close this introduction, let us mention that some first numerical investigations from the CRE-based strategy have been conducted for particular situations in [8,9], using input convex neural networks (ICNNs) to represent the learnt thermodynamic potentials (so-called NN-mCRE strategy). We also mention that the generalized standard material framework which is considered in this work is closely related to the (more general) GENERIC thermodynamic framework introduced in [27,67,70] and developed in solid mechanics in [63]; it was used in several works for data-driven modeling e.g., [14,32].
## Preliminaries: standard ill-posed problems and CRE concept

Material identification requires to solve problems with over-abundant data i.e., ill-posed computational problems. It is well-known that the classical formulation of such problems, without adding any regularization process, does not admit a solution. We will see that the Constitutive Relation Error (CRE) is able to give a solution, the constitutive relation being relaxed.

### Notations

The tested structure we shall study occupies the domain $\Omega \subset \mathbf{R}^3$ with boundary $\partial\Omega$, a typical point in $\Omega$ being denoted $\underline{M}$. The vectors, for instance displacements or velocities, are denoted by $\underline{U}$, $\underline{V}$, etc. Linear operators or second-order tensors, for instance strain or stress tensors, are written $\llap{\epsilon}$, $\sigma$, etc. More complex tensors or operators are written in bold text, such as the Hooke tensor which is written $\mathbf{K}$. This previous notation is also used for the different sets or spaces e.g., $\mathbf{R}^3$.

The Euclidean transpose is denoted by the superscript $^T$, so that the scalar product between vectors $\underline{V}$ and $\underline{W}$ reads $\underline{V}^T \underline{W}$. If there is no ambiguity, this scalar product can also be written $\underline{V} \cdot \underline{W}$. Work and power bring into play the trace operator such that:

$$ \mathrm{Tr}[\sigma \llap{\epsilon} \sigma] = \sum_{i,j \in \{1,2,3\}} [\sigma]_{ij} [\epsilon]_{ij} \quad (1) $$

In this paper, we shall consider small perturbations which means that the different configurations occupied by the structure can be replaced by the initial configuration $\Omega$. To fix ideas, the tested structure is placed in an environment characterized by:

* a displacement $\underline{U}_d$ on the part $\partial_1\Omega$ of the boundary $\partial\Omega$ ($\underline{U}_d \in [\mathbf{H}^{1/2}(\partial_1\Omega)]^3$);

* a surface force $\underline{F}_d$ on the part $\partial_2\Omega$ of the boundary $\partial\Omega$ ($\underline{F}_d \in [\mathbf{H}^{-1/2}(\partial_2\Omega)]^3$);

* a body force $\underline{f}_d$ over $\Omega$ (here we take $\underline{f}_d = 0$ for the sake of simplicity).

### Standard ill-posed problems

Commonly in material identification, the displacement field $\underline{U}_d$ is obtained from full-field measurements for which very efficient techniques are available today (see [33,71]). Moreover, boundary data is overabundant such that:

$$ \partial_1 \Omega \cup \partial_2 \Omega = \partial \Omega; \quad \partial_1 \Omega \cap \partial_2 \Omega = \partial_1 \Omega \neq \emptyset \quad (2) $$

which characterizes a standard ill-posed problem. These conditions are not mandatory for the use of the CRE concept, but they give certain properties that can help the calculation process. Of course, they can be extended to mixed boundary conditions. Incomplete conditions can also be introduced, but for this last case, the computational problem is no more standard.

In general, the part of the boundary where displacement full-field measurement takes place is also free, and therefore belongs to the part $\partial_{12}\Omega$ of the boundary where force is also known. A very common situation in identification problems is that the force is known over the entire boundary, eventually free on some parts, and the displacement field is measured on a part of the boundary. If the intersection of the parts of the boundary where the displacement is prescribed ($\partial_1\Omega$) and the force is given ($\partial_2\Omega$) is empty, i.e. $\partial_{12}\Omega = \emptyset$, the computational problem becomes well-posed if the boundary conditions cover the entire boundary; it has a solution for classical material models.
**Remark 1** We emphasize that a well-posed problem in elasticity, for example, is a problem for which the boundary conditions are such that it exists a unique solution for any constitutive relation. In other words, the knowledge of the boundary conditions does not give any information about the constitutive relation. This explains why one uses, in material identification, boundary conditions leading to an ill-posed computational problem.

Let us introduce $\mathcal{U}_{ad}^K$ the space of kinematically admissible displacement fields:

$$ \mathcal{U}_{ad}^K = \{ \underline{U} \text{ s.t. } \underline{U} \in [\mathbf{H}^1(\Omega)]^3, \underline{U}_{|\partial_1\Omega} = \underline{U}_d \} \quad (3) $$

This affine space can be written:

$$ \mathcal{U}_{ad}^K = \underline{U}_0^K + \mathcal{U}_{ad,0}^K \quad (4) $$

where $\underline{U}_0^K \in \mathcal{U}_{ad}^K$ is a particular admissible displacement field, and $\mathcal{U}_{ad,0}^K$ is the associated vectorial space.

Now, let us consider the space of statically admissible stress fields:

$$ \mathcal{S}_{ad}^S = \{ \sigma \text{ s.t. } \sigma \in [\mathbf{H}(\mathrm{div}, \Omega)]^3, \underline{\mathrm{div}}\sigma = 0 \text{ over } \Omega, \sigma \underline{n}_{|\partial_2\Omega} = \underline{F}_d \} \quad (5) $$

which can be defined from the principle of virtual works:

$$ \int_{\Omega} \mathrm{Tr}[\sigma \epsilon(\underline{U})] d\Omega = \int_{\partial_2 \Omega} \underline{F}_d \cdot \underline{U} dS \quad \forall \underline{U} \in \mathcal{U}_{ad,0}^S \quad (6) $$

thanks to the virtual displacement space:

$$ \mathcal{U}_{ad,0}^S = \{ \underline{U} \text{ s.t. } \underline{U} \in [\mathbf{H}^1(\Omega)]^3, \underline{U}_{|\partial_1\Omega \setminus \partial_{12}\Omega} = 0 \} \quad (7) $$

It follows the fundamental property of standard ill-posed problems:

**Property 1** For standard ill-posed problems, one has:

$$ \mathcal{U}_{ad,0}^K \subset \mathcal{U}_{ad,0}^S \quad (8) $$

## The Constitutive Relation Error

Let us consider the simple situation of elastic materials which may be linear or nonlinear; one has:

$$ \sigma = \mathbf{K}(\epsilon) \text{ over } \Omega \quad (9) $$

where the Hooke operator $\mathbf{K}$ is supposed to be known for the moment. An alternative formulation following J.J. Moreau's work [64] is to describe the elastic material behavior thanks to two potentials $\Psi(\epsilon)$ and $\Psi^*(\sigma)$ which are dual convex functions. This duality is defined by the work bilinear form:

$$ (\epsilon, \sigma) \in \mathbf{E} \times \mathbf{F} \mapsto \int_{\Omega} \mathrm{Tr}[\sigma \epsilon] d\Omega \in \mathbf{R} \quad (10) $$

One has the following well-known property:

**Property 2** Let $\Psi$ and $\Psi^*$ be two dual convex functions related to the work bilinear form. The CRE functional is then defined as:

$$ E_{CRE}^2(\underline{U}, \sigma) = \Psi(\epsilon(\underline{U})) + \Psi^*(\sigma) - \int_{\Omega} \mathrm{Tr}[\sigma \epsilon(\underline{U})] d\Omega \geq 0 \quad \forall (\epsilon, \sigma) \in \mathbf{E} \times \mathbf{F} \quad (11) $$

Moreover, (i) and (ii) are equivalent:

$$ (i) \quad \sigma = \mathbf{K}(\epsilon) \text{ over } \Omega \quad (12) $$

$$ (ii) \quad E_{CRE}(\underline{\nu}, \sigma) = 0 $$
$E_{CRE}$ is the global Constitutive Relation Error.

The classical formulation of the computational problem does not admit a solution for ill-posed problems. The main interest of the CRE method is to give a solution to such ill-posed problems by relaxing the constitutive relation. This is defined by:

$$ (\underline{\underline{U}}^K, \underline{\underline{\sigma}}^S) \in \mathcal{U}_{ad}^K \times \mathcal{S}_{ad}^S \equiv \underset{(\tilde{\underline{U}}^K, \tilde{\underline{\sigma}}^S) \in \mathcal{U}_{ad}^K \times \mathcal{S}_{ad}^S}{\arg \min} E_{CRE}^2(\tilde{\underline{U}}^K, \tilde{\underline{\sigma}}^S) \quad (13) $$

The value of the minimum $E_{CRE}^2(\underline{\underline{U}}^K, \underline{\underline{\sigma}}^S)$ is not equal to 0 in general for ill-posed computational problems.

Let us note that in the specific case of a linear elastic material behavior with $\underline{\underline{\sigma}} = \mathbf{K} \epsilon, \mathbf{K}$ being the Hooke elasticity tensor, corresponding potentials read $\Psi(\epsilon) = \frac{1}{2} \int_{\Omega} \mathrm{Tr}[\mathbf{K} \epsilon \epsilon] d\Omega$ and $\Psi^*(\underline{\underline{\sigma}}) = \frac{1}{2} \int_{\Omega} \mathrm{Tr}[\mathbf{K}^{-1} \underline{\underline{\sigma}} \underline{\underline{\sigma}}] d\Omega$, which yields:

$$ E_{CRE}^2(\underline{\underline{U}}, \underline{\underline{\sigma}}) = \frac{1}{2} \int_{\Omega} \mathrm{Tr}[(\underline{\underline{\sigma}} - \mathbf{K} \epsilon(\underline{\underline{U}}))\mathbf{K}^{-1}(\underline{\underline{\sigma}} - \mathbf{K} \epsilon(\underline{\underline{U}}))] d\Omega \quad (14) $$

### Computation of the optimal solution $(\tilde{\underline{U}}^K, \tilde{\underline{\sigma}}^S)$

In this section, we deal with how to practically compute the optimal solution of (13) for ill-posed computational problems.

From Property 1, we get the following decoupling property for any admissible pair $(\tilde{\underline{U}}^K, \tilde{\underline{\sigma}}^S) \in \mathcal{U}_{ad}^K \times \mathcal{S}_{ad}^S$:

$$ E_{CRE}^2(\tilde{\underline{U}}^K, \tilde{\underline{\sigma}}^S) = J_1(\tilde{\underline{U}}^K) + J_2(\tilde{\underline{\sigma}}^S) \quad (15) $$

with:

$$ J_1(\tilde{\underline{U}}^K) = \frac{1}{2} \int_{\Omega} \mathrm{Tr}[\epsilon(\tilde{\underline{U}}^K) \mathbf{K} \epsilon(\tilde{\underline{U}}^K)] d\Omega - \int_{\partial\Omega \setminus \partial_1 \Omega} \underline{\underline{E}}_d \cdot \tilde{\underline{U}}^K dS \\ + \frac{1}{2} \int_{\partial\Omega \setminus \partial_1 \Omega} \underline{\underline{E}}_d \cdot \underline{\underline{U}}_0^K dS \quad (16) $$

$$ J_2(\tilde{\underline{\sigma}}^S) = \frac{1}{2} \int_{\Omega} \mathrm{Tr}[\tilde{\underline{\sigma}}^S \mathbf{K}^{-1} \tilde{\underline{\sigma}}^S] d\Omega - \int_{\Omega} \mathrm{Tr}[\tilde{\underline{\sigma}}^S \epsilon(\underline{\underline{U}}_0^K)] d\Omega \\ + \frac{1}{2} \int_{\partial\Omega \setminus \partial_1 \Omega} \underline{\underline{E}}_d \cdot \underline{\underline{U}}_0^K dS $$

Indeed, for a given $\underline{\underline{U}}_0^K \in \mathcal{U}_{ad}^K$, we have $\tilde{\underline{U}}^K - \underline{\underline{U}}_0^K \in \mathcal{U}_{ad,0}^K \subset \mathcal{U}_{ad,0}^S$ which yields:

$$ \int_{\Omega} \mathrm{Tr}[\tilde{\underline{\sigma}}^S \epsilon(\tilde{\underline{U}}^K)] d\Omega = \int_{\Omega} \mathrm{Tr}[\tilde{\underline{\sigma}}^S \epsilon(\underline{\underline{U}}_0^K)] d\Omega + \int_{\partial\Omega \setminus \partial_1 \Omega} \underline{\underline{E}}_d \cdot (\tilde{\underline{U}}^K - \underline{\underline{U}}_0^K) dS \quad (17) $$

The decoupling property (15) is classical for well-posed computational problems (defining potential and complementary energies), but we just showed that it is also valid for standard ill-posed problems. It yields that the minimization (13) results in the solution of two independent minimization problems. Moreover, instead of solving the $J_2$-problem, we can compute its displacement formulation obtained thanks to a new dualization; defining

$$ J_1^S(\tilde{\underline{U}}^S) = \Psi(\epsilon(\tilde{\underline{U}}^S)) - \int_{\partial_2 \Omega} \underline{\underline{E}}_d \cdot \tilde{\underline{U}}^S d\Omega \quad (18) $$

and denoting by $\underline{\underline{U}}^S$ the minimizer of $J_1^S$ over $\mathcal{U}_{ad,0}^S$, we get $\underline{\underline{\sigma}}^S = \mathbf{K}(\epsilon(\underline{\underline{U}}^S))$.

Consequently, to compute the minimizer $(\underline{\underline{U}}^K, \underline{\underline{\sigma}}^S)$ of the CRE in $\mathcal{U}_{ad}^K \times \mathcal{S}_{ad}^S$, one computes the displacement pair $(\underline{\underline{U}}^K, \underline{\underline{U}}^S)$ which is easy. Of course, the classical finite element method can be used to perform these minimizations.
## Data-driven material modeling—elastic materials

The CRE-method is now developed for stable elastic materials in order to learn (9) (K is not given any more). First, we clearly specify the information which is available.

### Available information

#### Experimental data

We consider a series of tested structures $i \in \mathbf{N}_{test}$ which occupy the same domain $\Omega$. For each tested structure, the experimental data is:

$$ \underline{U}_d(i) \text{ over } \partial_1 \Omega; \quad \underline{F}_d(i) \text{ over } \partial_2 \Omega \qquad (19) $$

#### Additional knowledge from materials science

The considered environment is such that the studied material can be supposed stable. It follows that the behavior of the material can be described by two dual convex functions $\Psi$ and $\Psi^*$, the duality being defined by the work bilinear form. These two functions can be written:

$$ \Psi(\epsilon) = \int_{\Omega} \psi(\epsilon)d\Omega; \quad \Psi^*(\sigma) = \int_{\Omega} \psi^*(\sigma)d\Omega \qquad (20) $$

where $\psi$ and $\psi^*$ are two dual convex functions related to the local work bilinear form. These define the elastic material behavior as:

$$ \sigma = \frac{\partial \psi}{\partial \epsilon}; \quad \epsilon = \frac{\partial \psi^*}{\partial \sigma} \qquad (21) $$

Moreover, $\psi(\epsilon)$ is an energy and thus satisfies to:

$$ \psi \ge 0; \quad \psi(0) = 0; \quad \frac{\partial \psi}{\partial \epsilon}(0) = 0 \qquad (22) $$

The last relation means that the stress $\sigma$ is zero when the strain $\epsilon$ is also zero.

The corresponding admissible space for $\psi$ is denoted by $\psi$. It follows that the additional knowledge coming from materials science is reduced to:

$$ \psi \in \psi \qquad (23) $$

where $\psi$ is a convex subset.

### Optimal data-driven material model

#### Definition of the computed data-driven model

The computed data-driven model is the minimizer of the CRE for the family of tested structures. First, one defines the following CRE for a given material model (described by the thermodynamic potential $\psi \in \psi$):

$$ \bar{E}_{CRE}^2(\psi) = \frac{1}{N_{test}} \sum_{i \in \mathbf{N}_{test}} E_{CRE}^2(\psi; \underline{U}^K(i), \sigma^S(i)) \qquad (24) $$

with

$$ (\underline{U}^K(i), \sigma^S(i)) = \underset{(\bar{U}^K(i), \tilde{\sigma}^S(i)) \in \mathcal{U}_{ad}^K(i) \times \mathcal{S}_{ad}^S(i)}{\arg\min} E_{CRE}^2(\psi; \bar{U}^K(i), \tilde{\sigma}^S(i)) \qquad (25) $$

$\mathcal{U}_{ad}^K(i)$ and $\mathcal{S}_{ad}^S(i)$ being the admissible spaces related to problem $i \in \mathbf{N}_{test}$.

The computed data-driven model is thus determined by:

$$ \bar{\psi} \equiv \arg \min_{\psi \in \psi} \bar{E}_{CRE}^2(\psi) \qquad (26) $$
The final value $\bar{E}_{CRE}^2(\bar{\psi})$ reached by the CRE (averaged over the tested structures) at the end of the optimization process characterizes the quality of the optimal data-driven material model. It is different from 0 in general; in the very special case where this error is zero, it means the computed material model is compatible with the experimental data.

The validation of the calculated data-driven model, although a fundamental question, is beyond the scope of this paper. Nevertheless, a first idea is to describe the validation domain $\bar{s}$ from the set $s_{exp}$ of strain-stress pairs associated to the solutions of the i-problems with $\bar{\psi} \in \psi$. One has:

$$ \bar{s} = \left\{ (\epsilon, \sigma) \in \mathbf{R}^6 \times \mathbf{R}^6 \mid e_{CRE}(\bar{\psi}; \epsilon, \bar{\sigma}) \le \bar{\epsilon}, e_{CRE}(\bar{\psi}; \epsilon, \bar{\sigma}) \le \bar{\epsilon}, (\epsilon, \bar{\sigma}) \in s_{exp} \right\} \quad (27) $$

where the tolerance $\bar{\epsilon}$ satisfies $\bar{\epsilon} \ge \bar{E}_{CRE}(\bar{\psi})/|\Omega|$, and the local CRE is defined as:

$$ e_{CRE}^2(\bar{\psi}; \epsilon, \sigma) = \bar{\psi}(\epsilon) + \bar{\psi}^*(\sigma) - \mathrm{Tr}[\sigma\epsilon] \quad (28) $$

The experimental “points” should be sufficiently numerous and scattered so that the validation domain $\bar{\epsilon}$ contains a large ball of $\mathbf{R}^6$ centered at 0.

### Computation of the data-driven model

The problem to get $\bar{\psi} \in \psi$ being nonlinear, the computation method is iterative. At iteration $n + 1$, there are two stages:

1. K-S-stage at iteration $n + 1$ From the previous iterations, the following quantities have been computed: (i) $\psi_n \in \psi$; (ii) optimal fields $(\underline{U}_n^K(i), \sigma_n^S(i)) \in \mathcal{U}_{ad}^K(i) \times \mathcal{S}_{ad}^S(i)$, $i \in \mathbf{N}_{test}$, computed from $\psi_{n-1} \in \psi$. The K-S-stage consists in the computation of the solution of each i-problem, $i \in \mathbf{N}_{test}$, using the constitutive relation related to the convex function $\psi_n$. As we have seen in the previous section, this stage does not involve any serious difficulty for standard ill-posed computational problems; the classical finite element method can be used. In fact, we minimize the CRE for given $\psi_n$ to get optimal admissible fields:

$$ (\underline{U}_{n+1}^K(i), \sigma_{n+1}^S(i)) = \underset{(\underline{U}^K(i), \bar{\sigma}^S(i)) \in \mathcal{U}_{ad}^K(i) \times \mathcal{S}_{ad}^S(i)}{\arg \min} E_{CRE}^2(\psi_n; \underline{U}^K(i), \bar{\sigma}^S(i)) \quad (29) $$

and one has:

$$ \frac{1}{N_{test}} \sum_{i \in \mathbf{N}_{test}} E_{CRE}^2 (\psi_n; \underline{U}_{n+1}^K(i), \sigma_{n+1}^S(i)) \le \frac{1}{N_{test}} \sum_{i \in \mathbf{N}_{test}} E_{CRE}^2 ((\psi_n; \underline{U}_n^K(i), \sigma_n^S(i)) \quad (30) $$

The first term denoted $\bar{E}_{CRE}^2(\psi_n)$ can be interpreted as the CRE associated to the material model defined by $\psi_n \in \psi$. One has:

$$ \bar{E}_{CRE}^2(\psi_n) \le \bar{E}_{CRE}^2(\psi_{n-1}) \quad (31) $$

which indicates that the CRE decreases with the number of iterations. The iterative process is stopped when a stationary point is approximately reached. The final error value $\bar{E}_{CRE}(\bar{\psi})$ characterizes the quality of the computed data-driven material model.

2. $\psi$-stage at iteration $n + 1$ From the previous stage, one can define:

$$ \psi \in \psi \mapsto \bar{E}_{n+1}^2(\psi) = \frac{1}{N_{test}} \sum_{i \in \mathbf{N}_{test}} E_{CRE}^2 ((\psi; \underline{U}_{n+1}^K(i), \sigma_{n+1}^S(i)) \quad (32) $$
and $\psi_{n+1} \in \psi$ is thus computed as:

$$ \psi_{n+1} = \arg \min_{\psi \in \Psi} \bar{E}_{n+1}^2(\psi) \quad (33) $$

This corresponds here to a convex regression problem that is to find a convex function belonging to $\psi$ which fits the best possible the data, question which is detailed in the next section.

## Back to the $\psi$-stage at iteration $n + 1$

Convex regression is an old question with still a lot of research works nowadays. A first set of works refers to pioneer papers [34,69] where the convexity is obtained thanks to inequality constraints which are numerous and thus difficult to take into account. The detailed theory can be found in [51]. Among the numerous works in this direction, one can mention [4, 5, 65] and their references. Another set of works introduces the convex function being sought a priori, as a combination of elementary convex functions. Among these works, one can mention [30,61,74]. Input Convex Neural Networks (ICNN), introduced in [1] and widely used since then, follow the same representation of convex functions. We propose here to follow [74] where the convex function being sought is described as the envelop of quasi-affine convex functions. Using the Voigt notation where strain and stress tensors ($\epsilon, \sigma$) are represented by vectors ($\epsilon, \sigma$) of $\mathbf{R}^6$, the sought potential $\psi$ is represented as:

$$ \psi(\epsilon) = \max_{j \in \mathcal{J}} \psi_j(\epsilon); \quad \psi_j(\epsilon) = a_j + \underline{A}_j \cdot (\epsilon - \epsilon_j) + \frac{\tau}{2}(\epsilon - \epsilon_j)^T \mathbb{K} (\epsilon - \epsilon_j) \quad (34) $$

The positive symmetric matrix $\mathbb{K}$ and the curvature $\tau$ are given. For $j \in \mathcal{J}$, ($a_j, \underline{A}_j$) and $\epsilon_j$ are parameters. $\epsilon_j$ are part of the discretization points of the strain space.

The advantage of the representation (34), apart from the fact that it can be visualized, is that its Legendre-Fenchel transform can be defined explicitly; one has:

$$ \psi^*(\sigma) = \sup_{\epsilon \in \mathbb{R}^6} \left[ \underline{\sigma}^T \epsilon - \psi(\epsilon) \right] = \min_{j \in \mathcal{J}} \psi_j^*(\sigma) \quad (35) $$

with

$$ \psi_j^*(\sigma) = -a_j + \underline{\sigma}^T \epsilon_j + \frac{1}{2\tau} (\underline{\sigma} - \underline{A}_j)^T \mathbb{K}^{-1} (\underline{\sigma} - \underline{A}_j) \quad (36) $$

To belong to $\psi$, $\psi$ should in particular satisfy to:

$$ \psi(\underline{0}) = 0; \quad \frac{\partial \psi}{\partial \underline{\epsilon}}(\underline{0}) = \underline{0} \quad (37) $$

The last relation being equivalent to $\psi(\underline{0}) + \psi^*(\underline{0}) = 0$, one has:

$$ \begin{gathered} \max_{j \in \mathcal{J}} \left[ a_j - \underline{A}_j \cdot \epsilon_j + \frac{\tau}{2} \epsilon_j^T \mathbb{K} \epsilon_j \right] = 0 \quad (38) \\ \min_{j \in \mathcal{J}} \left[ -a_j + \frac{1}{2\tau} \underline{A}_j^T \mathbb{K}^{-1} \underline{A}_j \right] = 0 \end{gathered} $$

One also remarks that under the constraints (37):

$$ \forall \underline{\epsilon} \in \mathbb{R}^6, \quad \psi(\underline{\epsilon}) + \psi^*(\underline{0}) - 0 = \psi(\underline{\epsilon}) \ge 0 \quad (39) $$

and consequently (38) are the only constraints defining $\psi$.

It follows that the function $\bar{E}_{n+1}(\psi)$ introduced in (32) is explicitly defined in terms of the parameters ($a_j, \underline{A}_j, \underline{\epsilon}_j$) and consequently the problem to solve is a classical regression problem for which numerous numerical techniques are available. Of course, the optimal data-driven material model may be non-unique; consequently, a regularization (Lasso or Ridge) is welcome especially if the volume of data is relatively small.
## Particular case of linear materials

For linear elastic materials, it is remarkable that the minimization condition of the CRE functional can be written explicitly. Using the Voigt notation, the energy reads:

$$ \psi(\epsilon) = \frac{1}{2} \epsilon^T K \epsilon \qquad (40) $$

where

$$ K \in \mathcal{K} \equiv \{ K \text{ matrix } (6, 6) \text{ s.t. } K = K^T, K > 0 \} \qquad (41) $$

Here, the material model is thus defined by the Hooke tensor, the goal being to recover the optimal one.

The optimal $K$ is the minimizer over $\mathcal{K}$ of the CRE for the family of tested structures; one has following (24) and (26):

$$ \bar{E}_{CRE}^2(K) = \frac{1}{N_{test}} \sum_{i \in N_{test}} \min_{(\tilde{U}^K(i), \tilde{\sigma}^S(i)) \in \mathcal{U}_{ad}^K(i) \times \mathcal{S}_{ad}^S(i)} E_{CRE}^2(K; \tilde{U}^K(i), \tilde{\sigma}^S(i)) \qquad (42) $$

Denoting by $\epsilon^K$ and $\sigma^S$ the vectors corresponding respectively (in the Voigt notation) to optimal admissible fields $\epsilon(\tilde{U}^K)$ and $\sigma^S$, computed for given $K$ and therefore implicitly depending on $K$, one has:

$$ \bar{E}_{CRE}^2(K) = \int_{\Omega} \left( \frac{1}{2} \mathrm{Tr}[KQ^K] + \frac{1}{2} \mathrm{Tr}[K^{-1}Q^S] - \frac{1}{N_{test}} \sum_{i \in N_{test}} \underline{\sigma}^S(i)^T \underline{\epsilon}^K(i) \right) d\Omega \qquad (43) $$

with

$$ Q^K = \frac{1}{N_{test}} \sum_{i \in N_{test}} \epsilon^K(i) \epsilon^K(i)^T; \quad Q^S = \frac{1}{N_{test}} \sum_{i \in N_{test}} \underline{\sigma}^S(i) \underline{\sigma}^S(i)^T \qquad (44) $$

The differential of $\bar{E}_{CRE}$ should be equal to zero, which yields (using $\delta K = -K \delta K^{-1} K$):

$$ 0 = -\frac{1}{2}\mathrm{Tr}[\delta K^{-1}KQ^KK] + \frac{1}{2}\mathrm{Tr}[\delta K^{-1}Q^S] + \left[ \frac{1}{2}\mathrm{Tr}[K\delta Q^K] + \frac{1}{2}\mathrm{Tr}[K^{-1}\delta Q^S] - \delta \frac{1}{N_{test}} \sum_{i \in N_{test}} \underline{\sigma}^S(i)\underline{\epsilon}^K(i)^T \right] \qquad (45) $$

One can note that the last term between brackets is equal to zero because $(\epsilon^K(i), \sigma^S(i))$ is the minimizer of the CRE for the $i$-problem. Finally, one gets:

**Property 3** For linear elastic materials, the minimization condition of the CRE is:

$$ KQ^K = Q^S; \quad K \in \mathcal{K} \qquad (46) $$

Let us come back to the $\psi$-stage at iteration $n+1$. One has:

$$ KQ_{n+1}^K = Q_{n+1}^S; \quad K \in \mathcal{K} \qquad (47) $$

so that

$$ \left( \left[ Q_{n+1}^K \right]^{1/2} K \left[ Q_{n+1}^K \right]^{1/2} \right) \left( \left[ Q_{n+1}^K \right]^{1/2} K \left[ Q_{n+1}^K \right]^{1/2} \right) = \left[ Q_{n+1}^K \right]^{1/2} Q_{n+1}^S \left[ Q_{n+1}^K \right]^{1/2} \qquad (48) $$

Consequently, if $Q_{n+1}^K$ and $Q_{n+1}^S$ are strictly positive:

$$ K^{n+1} = \left[ Q_{n+1}^K \right]^{-1/2} \left( \left[ Q_{n+1}^K \right]^{1/2} Q_{n+1}^S \left[ Q_{n+1}^K \right]^{1/2} \right)^{1/2} \left[ Q_{n+1}^K \right]^{-1/2} \qquad (49) $$
The case where matrices $Q_{n+1}^K$ and $Q_{n+1}^S$ are singular i.e., when the strain and stress spaces are not fully explored by the tests, is more complicated. We will omit the index $n + 1$. Let $\pi_K$ and $\pi_S$ be the orthogonal projectors on the kernels of $Q^K$ and $Q^S$, respectively. To define the part of $K$ which is not defined by the minimization condition, we use a regularization associated to the reference Hooke tensor $K_0$ (strictly positive).

One remarks first from (48) that:

$$ \|Q^{K^{1/2}} \mathbb{K} Q^{K^{1/2}} \underline{W}\|^2 = \underline{W}^T Q^{K^{1/2}} Q^S Q^{K^{1/2}} \underline{W} \quad (50) $$

If $(\mathbb{I} - \pi_K) \underline{W} \in \text{Ker}(Q^S)$, the part $(\mathbb{I} - \pi_K) \mathbb{K}(\mathbb{I} - \pi_K)$ of the operator $\mathbb{K}$ is singular; one has:

$$ (\mathbb{I} - \pi_K) \mathbb{K}(\mathbb{I} - \pi_K) = Q^{K^{-1/2}} \left[ Q^{K^{1/2}} Q^S Q^{K^{1/2}} \right]^{1/2} Q^{K^{-1/2}} \quad (51) $$

To get a regular operator, let us consider the orthogonal projection $(\mathbb{I} - \bar{\pi}_K)$ on the subspace $\text{Ker}^\perp(Q^K) \cap \text{Ker}(Q^S)$; one takes to get a non-singular $\mathbb{K}$:

$$ (\mathbb{I} - \pi_K) \mathbb{K}(\mathbb{I} - \pi_K) = Q^{K^{-1/2}} \left[ Q^{K^{1/2}} Q^S Q^{K^{1/2}} \right]^{1/2} Q^{K^{-1/2}} + (\mathbb{I} - \bar{\pi}_K) K_0 (\mathbb{I} - \bar{\pi}_K) \quad (52) $$

Moreover, one has:

$$ (\mathbb{I} - \pi_K) \mathbb{K} \pi_K = Q^{K^{-1}} \left[ (\mathbb{I} - \pi_K) \mathbb{K}(\mathbb{I} - \pi_K) \right]^{-1} Q^S \pi_K \quad (53) $$

To end, we take:

$$ \pi_K \mathbb{K} \pi_K = \pi_K K_0 \pi_K \quad (54) $$

## Comparison with the pioneer data-driven approach

Let us consider the very particular situation where strain and stress are uniform for the family $\mathbf{N}_{test}$ of tested structures. It follows that one can easily define the Experimental Constitutive Manifold $\Gamma_{exp}$:

$$ \Gamma_{exp} = \left\{ (\epsilon^K(i), \omega^S(i)) \in \mathbf{R}^6 \times \mathbf{R}^6, i \in \mathbf{N}_{test} \right\} \quad (55) $$

The pioneer work [41] starts with the Experimental Constitutive Manifold which is in general not explicitly defined in terms of experimental data. To get the data-driven material model, an interpolation operator is added, for instance using a nearest neighbors (k-NN) algorithm. The present work is different in the sense that from materials science, the available mathematical shapes of the Experimental Constitutive Manifold are given. For linear elastic materials, $\Gamma_{exp}$ should be close to a linear manifold that is computed thanks to the CRE. One gets for tests covering the strain and stress spaces:

$$ \mathbb{K} = [Q^K]^{-1/2} \left( [Q^K]^{1/2} Q^S [Q^K]^{1/2} \right)^{1/2} [Q^K]^{-1/2} \quad (56) $$

with

$$ Q^K = \frac{1}{N_{test}} \sum_{i \in \mathbf{N}_{test}} \epsilon^K(i) \epsilon^K(i)^T; \quad Q^S = \frac{1}{N_{test}} \sum_{i \in \mathbf{N}_{test}} \sigma^S(i) \sigma^S(i)^T \quad (57) $$

In the case of nonlinear but stable elastic materials, $\Gamma_{exp}$ is nonlinear and should be close to a manifold defined by the convex potential $\psi \in \psi$. On has by minimizing the CRE:

$$ \bar{\psi} = \arg \inf_{\psi \in \psi} \frac{1}{N_{test}} \sum_{i \in \mathbf{N}_{test}} \left[ \psi(\epsilon^K(i)) + \psi^*(\omega^S(i)) - \text{Tr}[\omega^S(i)\epsilon^K(i)] \right] \quad (58) $$

## Data-driven material modeling—stable elasto-(visco-)plastic materials

The CRE-method is now set up for the data-driven modeling of stable elasto-(visco-)plastic materials. There are two key questions for which a solution is given:
* the mathematical shape of the most general data-driven model compatible with knowledge from physics and materials science;

* the computation of standard ill-posed problems for elasto-(visco-)plastic materials.

## Experimental data

We use similar notations as in the previous sections, but introducing now the time domain. We have a series of tested structures $i \in \mathbf{N}_{test}$ over the time interval $[0, T]$. Data obtained from measurements is:

$$ \underline{U}_d(i) \text{ over } \partial_1 \Omega \times [0, T]; \quad \underline{E}_d(i) \text{ over } \partial_2 \Omega \times [0, T] \qquad (59) $$

The i-problems ($i \in \mathbf{N}_{test}$) are standard ill-posed problems for which the admissible spaces for the displacement-stress pair $(\underline{U}^K, \underline{\omega}^S)$ are denoted $\mathcal{U}_{ad}^{K[0,T]}$ and $\mathcal{S}_{ad}^{S[0,T]}$.

## Mathematical shape of data-driven models

The goal in this section is to derive the mathematical shape of the most general available model, the elastic stiffness being known. We use here the internal variables framework which is very general and common in materials science; this is the basis of the material modeling developed here for stable elasto-(visco-)plastic materials.

### The different internal variables

The state of the material at time $t$ is defined by the value of the following set:

$$ (\epsilon_p, \mathbf{X}, \omega, \mathbf{Y}) \qquad (60) $$

where $\epsilon_p (= \epsilon - K^{-1}\omega)$ and $\omega$ are observable variables, while the pair $(\mathbf{X}, \mathbf{Y})$ gathers hidden internal variables. $\mathbf{Y}$ is the force associated to $\mathbf{X}$ and both belong to $\mathbf{R}^q$. In a data-driven mathematical model, $\mathbf{X}$ and $\mathbf{Y}$ are a priori unknown as well as their dimension $q$.

### The classical formulation

For the sake of simplicity, one assumes isothermal conditions. Following the first principle of thermodynamics of irreversible processes, we consider a Helmholtz free energy $e$ which depends on the state of the material defined by the elastic strain $\epsilon_e$, the hidden variable $\mathbf{X} \in \mathbf{R}^q$ and the cumulated inelastic strain $p$:

$$ e(\epsilon_e, \mathbf{X}, p) \qquad (61) $$

$\mathbf{X}$ has often a precise interpretation; for example, in the case of kinematic hardening, $\mathbf{X}$ is associated in the classical models to the center of the elasticity domain.

As classically done in materials science, one introduces the following decoupling hypothesis ($H_1$) used in this section:

$$ e = \frac{1}{2} \mathrm{Tr}[\mathbf{K}\epsilon_e\epsilon_e] + \psi(\mathbf{X}, p) \qquad (62) $$

It follows the state equations defining the forces associated to $(\epsilon_e, \mathbf{X}, p)$:

$$ \omega = \mathbf{K}\epsilon_e; \quad \mathbf{Y} = \frac{\partial\psi}{\partial\mathbf{X}}; \quad R = \frac{\partial\psi}{\partial p} \qquad (63) $$

The dissipation is equal to:

$$ \begin{aligned} \omega &= \mathrm{Tr}[\omega\dot{\epsilon}] - \dot{e} \\ &= \mathrm{Tr}[\omega\dot{\epsilon}_p] - \dot{\psi} \\ &= \mathrm{Tr}[\omega\dot{\epsilon}_p] - \mathbf{Y} \cdot \dot{\mathbf{X}} - R\dot{p} \end{aligned} \qquad (64) $$
The second principle of thermodynamics of irreversible processes indicates that the dissipation $\omega$ is positive or null for all admissible histories of the material. To satisfy this constraint, one adds the state evolution laws:

$$ \begin{bmatrix} \dot{\epsilon}_p \\ -\dot{\mathbf{X}} \\ -\dot{p} \end{bmatrix} = \mathbf{B} \begin{pmatrix} \omega \\ \mathbf{Y} \\ R \end{pmatrix} \qquad (65) $$

with $\epsilon_p = 0$, $\mathbf{X} = 0$ and $p = 0$ at $t = 0$, and $\mathbf{B}$ is positive.

### The normal formulation

The hidden internal variables $(\mathbf{X}, \mathbf{Y}) \in \mathbf{R}^q \times \mathbf{R}^q$ are not intrinsic in the sense that another pair $(\mathbf{X}', \mathbf{Y}') \in \mathbf{R}^q \times \mathbf{R}^q$ can lead to the same observable quantities i.e. the same material model. This legitimates (as done in [47]) a choice more remarkable than the others, the so-called normal formulation. It results from a transformation of hidden variables, from $(\mathbf{X}, \mathbf{Y})$ to $(\tilde{\mathbf{X}}, \tilde{\mathbf{Y}})$, which leads to an equivalence between the kinematic and static variables, of the form:

$$ \tilde{\mathbf{X}} = \Lambda \tilde{\mathbf{Y}} \qquad (66) $$

where $\Lambda$ is a positive symmetric constant linear operator. Here, we take $\Lambda = I$. This formulation is not unique, and the following property introduces a particular one where the only assumption is the decoupling hypothesis $(H_1)$. Moreover, this normal formulation preserves the energy and the dissipation.

**Property 4** Let the Helmholtz free energy be:

$$ e(\epsilon_e, \mathbf{X}, p) = \frac{1}{2} \mathrm{Tr}[\mathbf{K}\epsilon_e\epsilon_e] + \psi(\mathbf{X}, p) \qquad (67) $$

where $\psi$ is a convex positive function, twice differentiable such that:

$$ \psi(\mathbf{0}, 0) = 0; \quad \frac{\partial \psi}{\partial \mathbf{X}}(\mathbf{0}, 0) = 0; \quad \frac{\partial \psi}{\partial p}(\mathbf{0}, 0) = 0 \qquad (68) $$

There is an internal variables transform:

$$ \tilde{\mathbf{X}} = \mathbf{Q}_v(\mathbf{X}, p); \quad \tilde{p} = q_v(\mathbf{X}, p) \qquad (69) $$

such that the energy can be written with the new variables $(\tilde{\mathbf{X}}, \tilde{p})$:

$$ e = \frac{1}{2} \mathrm{Tr}[\mathbf{K}\epsilon_e\epsilon_e] + \frac{1}{2}\tilde{\mathbf{X}} \cdot \tilde{\mathbf{X}} + \frac{1}{2}\tilde{p}^2 \qquad (70) $$

The dissipation is also preserved. The force $\tilde{\mathbf{Y}}$ associated to $\tilde{\mathbf{X}}$ is such that $\tilde{\mathbf{Y}} = \tilde{\mathbf{X}}$ and $\tilde{p}$ can be taken in $\mathbf{R}^+$.

**Proof** Let $\mathbf{Z} = \begin{bmatrix} \mathbf{X} \\ p \end{bmatrix} \in \mathbf{R}^{q+1}$ and

$$ \bar{\psi}(\mathbf{Z}) = \mathbf{Z}^T \left[ \int_0^1 d\mu \int_0^1 \lambda d\lambda \frac{\partial^2 \psi}{\partial \mathbf{Z}^2}(\lambda \mu \mathbf{Z}) \right] \mathbf{Z} \qquad (71) $$

First, we prove that $\bar{\psi}(\mathbf{Z}) = \psi(\mathbf{Z})$. One has:

$$ \bar{\psi}(\mathbf{Z}) = \int_0^1 d\mu \frac{d}{d\mu} \Gamma(\mu) \qquad (72) $$
with $\Gamma(\mu) = \int_{0}^{1} d\lambda \frac{\partial \psi}{\partial Z}(\mu\lambda Z)Z$. It follows:

$$ \overline{\psi}(Z) = \Gamma(1) - \Gamma(0) = \int_{0}^{1} d\lambda \frac{\partial \psi}{\partial Z}(\lambda Z)Z \quad (73) $$

Finally, one has:

$$ \overline{\psi}(Z) = \psi(Z) - \psi(0) = \psi(Z) \quad (74) $$

Let us introduce the operator

$$ A = 2 \int_{0}^{1} \int_{0}^{1} \lambda d\lambda d\mu \frac{\partial^2 \psi}{\partial Z^2}(\lambda \mu Z) \quad (75) $$

It is symmetric and also positive, $\psi$ being a convex function. Consequently, one can define the operator $A^{1/2}$ and the new state variables are:

$$ \begin{bmatrix} \tilde{X} \\ \tilde{p} \end{bmatrix} = A^{1/2} \begin{bmatrix} X \\ p \end{bmatrix} \quad (76) $$

with $\psi(\tilde{X}, \tilde{p}) = \frac{1}{2}\tilde{X} \cdot \tilde{X} + \frac{1}{2}\tilde{p}^2$.

The energy being preserved, it is also true for the dissipation. The operator **B** which defines the state evolution laws is of course modified, the relation between the observable variables being preserved. □

The normal formulation has serious advantages. For stable materials with a convex energy, it allows to divide by two the number of hidden internal variables and to get an energy parametrized by the Hooke tensor **K** alone.

Elasto-viscoplasticity is characterized by the fact that the components of the generalized strain rate, i.e. $\epsilon_p, \dot{\tilde{X}}, \dot{\tilde{p}}$, are null at the same time. Consequently, the most general mathematical model can be written:

$$ \begin{align} \epsilon_p &= g(\sigma, \tilde{X}, \tilde{p})\dot{\tilde{p}} \notag \\ \dot{\tilde{X}} &= h(\sigma, \tilde{X}, \tilde{p})\dot{\tilde{p}} \notag \\ \dot{\tilde{p}} &= f(\sigma, \tilde{X}, \tilde{p}) \ge 0 \tag{77} \\ \epsilon_p = 0, \quad \tilde{X} = 0, \quad \tilde{p} = 0 \quad \text{at } t = 0 \notag \end{align} $$

This model is defined by three functions and the number $q$ of hidden variables. To fulfill the second principle of thermodynamics of irreversible processes, one should satisfy:

$$ \mathrm{Tr}[\sigma \epsilon_p] - \tilde{X} \cdot \dot{\tilde{X}} - \tilde{p} \cdot \dot{\tilde{p}} \ge 0 \quad (78) $$

For elastoplastic materials for which the behavior is independent of the loading rate, the third equation of (77) is replaced by:

$$ \dot{\tilde{p}} \bar{f}(\sigma, \tilde{X}, \tilde{p}) = 0; \quad \dot{\tilde{p}} \ge 0; \quad \bar{f} \le 0 \quad (79) $$

**Remark 2** We have shown in [50] that the functions **g**, **h**, **f** should be single-valued to get a consistent material model with (77) and (78). This univocity property is the basis of the construction method of the hidden internal state variables from data introduced in [24,50]. Its interest is its generality. However, we will see that the CRE enables to solve the same problem with a more robust approach but less general in the sense that a CRE must be defined beforehand, constraint that we are considering today as relatively weak.
### The subfamily of standard materials

The decoupling hypothesis ($H_1$) and the convexity property of the energy $e$ belong to commonly accepted assumptions in materials science as they have been validated by experiments for stable materials. Consequently, they will be systematically taken into account a priori in the following. Here, we go one step further considering stable materials which, for most of them, can be described by a standard material model. Such a family, very common in materials science, has been validated for numerous materials. It has been introduced by B. Halphen and Q.S. Nguyen [29].

Precisely, a standard material model is defined by two potentials $\phi^*(\sigma, \tilde{X}, \tilde{p})$ and $\phi(\epsilon_p, -\dot{\tilde{X}}, -\dot{\tilde{p}})$ which are two dual convex functions related to the dissipation bilinear form:

$$ ((\epsilon_p, -\dot{\tilde{X}}, -\dot{\tilde{p}}), (\sigma, \tilde{X}, \tilde{p})) \in \mathbf{R}^{6+q+1} \times \mathbf{R}^{6+q+1} \mapsto \mathrm{Tr}[\sigma \epsilon_p] - \tilde{X} \cdot \dot{\tilde{X}} - \tilde{p} \cdot \dot{\tilde{p}} \quad (80) $$

One should also satisfy the second principle of thermodynamics of irreversible processes, that corresponds to:

$$ \phi^*(0, \mathbf{0}, 0) = 0; \quad \phi^* \ge 0 \quad (81) $$

Such a material model depends only on one convex function and it is compatible with the two principles of thermodynamics. It also depends on the dimension $q$ of the vector of hidden internal variables.

The local CRE then reads:

$$ e_{CRE}^2 = \phi(\epsilon_p, -\dot{\tilde{X}}, -\dot{\tilde{p}}) + \phi^*(\sigma, \tilde{X}, \tilde{p}) - \mathrm{Tr}[\sigma \epsilon_p] + \tilde{X} \cdot \dot{\tilde{X}} + \tilde{p} \cdot \dot{\tilde{p}} \quad (82) $$

and the global CRE over the time-space domain $[0, T] \times \Omega$ is equal to:

$$ E_{CRE}^2 = \int_0^T dt \int_{\Omega} d\Omega \eta(t) e_{CRE}^2 \quad (83) $$

where the weight $\eta(t)$ may be chosen as 1 or $(1 - t/T)$.

Let us return to the normal formulation and the following question: is there a normal formulation for a standard initial model that retains the property of being standard? The answer has been given in [47]; one has:

**Property 5** A standard material model defined by the convex potential $\phi^*(\sigma, Z)$ and the inelastic part of the energy $\psi(Z)$ admits a normal standard formulation is there exists an operator $\mathbf{R}$ (defining a new set of hidden internal variables) such that:

$$ \tilde{Z} = \begin{bmatrix} \tilde{X} \\ \tilde{p} \end{bmatrix} = \mathbf{R}(Z) \qquad \left(\frac{\partial \tilde{Z}}{\partial Z}\right)^T \left(\frac{\partial \tilde{Z}}{\partial Z}\right) = \frac{\partial}{\partial Z} \left(\frac{\partial \psi}{\partial Z}\right)^T \quad (84) $$

$$ \phi^* \left( \sigma, \left( \frac{\partial \psi}{\partial Z} \right)^T (\mathbf{R}^{-1}(\tilde{Z})) \right) \text{ is a convex function} $$

This normal formulation is different from the one introduced in “The normal formulation” section. The conditions of Property 5 are much more common than they appear at first. They are satisfied by classical parametrized models of materials science [47]; an illustration is given a little further on.
One can remark that the main condition:

$$ \left(\frac{\partial \tilde{Z}}{\partial \tilde{Z}}\right)^T \left(\frac{\partial \tilde{Z}}{\partial \tilde{Z}}\right) = \frac{\partial}{\partial \tilde{Z}} \left(\frac{\partial \psi}{\partial \tilde{Z}}\right)^T \qquad (85) $$

involves a deformation. The right-hand side should be positive which implies the convexity of the energy. An approximation can be given if we neglect the square of the nonlinear part of $\mathbf{R}(\mathbf{Z})$ denoted $\mathbf{R}^*(\mathbf{Z})$; one has:

$$ \tilde{Z} = \mathbf{R}(\mathbf{Z}) = \mathbf{A}_0 \mathbf{Z} + \mathbf{R}^*(\mathbf{Z}) \qquad (86) $$

It follows with $\mathbf{\Delta} = \frac{\partial}{\partial \mathbf{Z}} \mathbf{R}^*(\mathbf{Z})$:

$$ \mathbf{A}_0^T \mathbf{\Delta} + \mathbf{\Delta}^T \mathbf{A}_0 + \mathbf{A}_0^T \mathbf{A}_0 = \frac{\partial}{\partial \mathbf{Z}} \left( \frac{\partial \psi}{\partial \mathbf{Z}} \right)^T \qquad (87) $$

One gets the operator $\mathbf{R}$ as:

$$ \tilde{Z} = \frac{1}{2} \left[ \left( \mathbf{A}_0^{-1} \right)^T \left( \frac{\partial \psi}{\partial \mathbf{Z}} \right)^T + \mathbf{A}_0 \mathbf{Z} \right] \qquad (88) $$

which shows that if we neglect the square of the nonlinear part of $\frac{\partial \psi}{\partial \mathbf{Z}})^T$, the material model is both normal and standard.

Finally, as a conclusion of the previous developments, we consider that the family of normal and standard models provides a reasonable framework for the data-driven modeling of stable elasto-(visco-)plastic materials.

### A classical model in materials science

Let us consider first the standard version of the Chaboche-Marquis elasto-(visco-)plastic model (see [47,52]). We will see that this model admits a version which is both standard and normal.

The material is supposed to be isotropic; the internal hidden variable is a second order tensor $\mathbf{X}$ with Tr[$\mathbf{X}$] = 0, its associated force being $\mathbf{Y}$ with Tr[$\mathbf{Y}$] = 0. The dissipation potential in viscoplasticity is with classical notations:

$$ \phi^*(\sigma_D, \mathbf{Y}, \mathbf{R}) = \frac{k}{n+1} \langle z \rangle_+^{n+1} \quad n > 0 \qquad (89) $$

with

$$ \begin{aligned} z &= |\sigma_D - \mathbf{Y}| + \frac{1}{2} a |\mathbf{Y}|^2 - \mathbf{R} \\ \mathbf{Y} &= c\mathbf{X} \\ R &= \gamma(p) \\ p &= \int_0^T |\dot{\epsilon}_p| dt \end{aligned} \qquad (90) $$

$k, n, a, c$ are positive material parameters. The threshold function:

$$ p \in \mathbf{R}^+ \mapsto \gamma(p) \in \mathbf{R}^+ \qquad (91) $$

is also a material function which is supposed to be concave and increasing.

To get a normal standard formulation, let us introduce:

$$ \begin{aligned} \tilde{\mathbf{X}} &= c^{1/2} \mathbf{X} \\ \tilde{p} &= \int_0^p dp \left[ \frac{\partial \gamma}{\partial p} \right]^{1/2} \\ R &= \ell(\tilde{p}) \end{aligned} \qquad (92) $$
One has $\tilde{\phi}^*(\sigma_D, \tilde{\mathbf{X}}, \tilde{p}) = \phi^*(\sigma_D, c^{1/2}\mathbf{X}, \ell(\tilde{p}))$ and one can check that $\tilde{\phi}^*$ is convex.

**Standard material model for elastoplastic materials**

Let us consider the case of elastoplastic materials. The potential $\phi^*$ is the indicator function of the elastic domain which is convex; one has:

$$ \phi^*(\sigma, \tilde{\mathbf{X}}, \tilde{p}) = \begin{cases} 0 & \text{if } \bar{f}(\sigma, \tilde{\mathbf{X}}, \tilde{p}) \le 0 \\ +\infty & \text{otherwise} \end{cases} \qquad (93) $$

with $\bar{f}(0, 0, 0) \le 0$.

This mathematical model depends also on one convex function $\bar{f}$. This function can be simplified because we know from materials science that the elasticity domain grows with $\tilde{p}$ which can be interpreted as a material time, i.e.:

$$ \bar{f}_{,\tilde{p}} < 0 \quad (H_2) \qquad (94) $$

It follows the property:

**Property 6** Under the $(H_2)$ hypothesis and if $\bar{f}$ is convex and piecewise regular, $\bar{f}$ can be replaced by:

$$ a(\sigma, \tilde{\mathbf{X}}) - \tilde{p} \qquad (95) $$

where a is a convex function such that $a(0, 0) \le 0$.

**Proof** From $\bar{f}(\sigma, \tilde{\mathbf{X}}, \tilde{p}) = 0$, one gets

$$ \tilde{p} = a(\sigma, \tilde{\mathbf{X}}) \qquad (96) $$

and we have:

$$ \bar{f}(\sigma, \tilde{\mathbf{X}}, a(\sigma, \tilde{\mathbf{X}})) = 0 \qquad (97) $$

The second differential of $\bar{f}(\sigma, \tilde{\mathbf{X}}, a(\sigma, \tilde{\mathbf{X}}))$ is:

$$ 0 = [\delta^2 \bar{f}](\sigma, \tilde{\mathbf{X}}, \tilde{p})_{\tilde{p}=a(\sigma, \tilde{\mathbf{X}})} + \frac{\partial \bar{f}}{\partial p}[\delta^2 a](\sigma, \tilde{\mathbf{X}}) \qquad (98) $$

It follows $\delta^2 a \ge 0$ and consequently the convexity of $a$. Moreover, from $\bar{f}(0, 0, 0) \le 0$, $p \ge 0$, one gets:

$$ a(0, 0) \le 0 \qquad (99) $$

Finally, the data-driven mathematical model compatible with the knowledge from physics and materials science is defined by the convex function $a$ alone, such that $a(0, 0) \le 0$. □

**Remark 3** Some more assumptions can be introduced for particular materials. For example, for metallic materials, the plastic incompressibility is commonly used. Isotropy or orthotropy properties can also be added and allow important complexity reduction.

It follows that the dual potential $\phi$ is:

$$ \begin{aligned} \phi(\epsilon_p, -\dot{\tilde{\mathbf{X}}}, -\dot{\tilde{p}}) &= \sup_{\sigma, \tilde{\mathbf{X}}, \tilde{p}, a(\sigma, \tilde{\mathbf{X}}) \le \tilde{p}} \left[ \mathrm{Tr}[\sigma \epsilon_p] - \tilde{\mathbf{X}} \cdot \dot{\tilde{\mathbf{X}}} - \tilde{p} \cdot \dot{\tilde{p}} \right] \\ &= \sup_{\sigma, \tilde{\mathbf{X}}, \tilde{p}, a(\sigma, \tilde{\mathbf{X}}) \le \tilde{p}} \left[ \dot{\tilde{p}}(\mathrm{Tr}[\sigma \epsilon_p / \dot{\tilde{p}}]) - \tilde{\mathbf{X}} \cdot \dot{\tilde{\mathbf{X}}} / \dot{\tilde{p}} - a(\sigma, \tilde{\mathbf{X}}) + \dot{\tilde{p}}(a - \tilde{p}) \right] (100) \\ &= \dot{\tilde{p}} a^*(\epsilon_p / \dot{\tilde{p}}, -\dot{\tilde{\mathbf{X}}} / \dot{\tilde{p}}) + 0 \end{aligned} $$
The local CRE can be written:

$$e_{CRE}^2 = \dot{p} \left[ a^*(\epsilon_p/\dot{p}, -\dot{\mathbf{X}}/\dot{p}) + a(\sigma, \dot{\mathbf{X}}) - \mathrm{Tr}[\sigma \epsilon_p/\dot{p}] + \dot{\mathbf{X}} \cdot \dot{\mathbf{X}}/\dot{p} \right] + \dot{p} \left[ \tilde{p} - a(\sigma, \dot{\mathbf{X}}) \right] \quad (101)$$

where the quantities between brackets are positive or null.

**Standard material model for elasto-viscoplastic materials**

A common assumption in materials science is to consider viscoplasticity as an extension or regularization of plasticity. The elastic domain being described by the convex function a, the potential is defined as:

$$\phi^*(\sigma, \tilde{\mathbf{X}}, \tilde{p}) = \gamma(\langle a(\sigma, \tilde{\mathbf{X}}) - \tilde{p} \rangle_+) \quad (102)$$

where $\gamma$ is an increasing convex function defined over $\mathbf{R}^+$ such that $\gamma(0) = 0$. The potential $\phi^*$ is thus convex, positive and satisfies:

$$\phi^*(0, \mathbf{0}, 0) = 0 \quad (103)$$

To get the dual potential $\phi$, one has to take the Legendre-Fenchel transform of a combination of convex functions named composite convex functions for which results are given in [36]. Here, one has:

$$\begin{aligned}
\phi(\epsilon_p, -\dot{\mathbf{X}}, -\dot{p}) &= \sup_{\sigma, \tilde{\mathbf{X}}, \tilde{p} \in \mathbf{R}^{6+q} \times \mathbf{R}^+} \left[ \mathrm{Tr}[\sigma \epsilon_p] - \tilde{\mathbf{X}} \cdot \dot{\mathbf{X}} - \tilde{p} \cdot \dot{p} - \gamma(\langle a(\sigma, \tilde{\mathbf{X}}) - \tilde{p} \rangle_+) \right] \\
&= \sup_{\sigma, \tilde{\mathbf{X}}, \tilde{p} \in \mathbf{R}^{6+q} \times \mathbf{R}^+} \left[ (\dot{p}(a(\sigma, \tilde{\mathbf{X}}) - \tilde{p}) - \gamma(\langle a(\sigma, \tilde{\mathbf{X}}) - \tilde{p} \rangle_+)) \right. \\
& \qquad \left. + \dot{p}(a^*(\epsilon_p/\dot{p}, -\dot{\mathbf{X}}/\dot{p}) + a(\sigma, \tilde{\mathbf{X}}) - \mathrm{Tr}[\sigma \epsilon_p/\dot{p}] + \tilde{\mathbf{X}} \cdot \dot{\mathbf{X}}/\dot{p}) \right] \\
&= \gamma^*(\dot{p}) + \dot{p} a^*(\epsilon_p/\dot{p}, -\dot{\mathbf{X}}/\dot{p})
\end{aligned} \quad (104)$$

It follows that the local CRE is equal to:

$$\begin{aligned}
e_{CRE}^2 = \dot{p} \left[ a(\sigma, \tilde{\mathbf{X}}) + a^*(\epsilon_p/\dot{p}, -\dot{\mathbf{X}}/\dot{p}) - \mathrm{Tr}[\sigma \epsilon_p/\dot{p}] + \tilde{\mathbf{X}} \cdot \dot{\mathbf{X}}/\dot{p} \right] \\
+ \left[ \gamma^*(\dot{p}) + \gamma(\langle a(\sigma, \tilde{\mathbf{X}}) - \tilde{p} \rangle_+) - \dot{p} \langle a(\sigma, \tilde{\mathbf{X}}) - \tilde{p} \rangle_+ \right] \\
+ \left[ \dot{p} \langle a(\sigma, \tilde{\mathbf{X}}) - \tilde{p} \rangle_- \right]
\end{aligned} \quad (105)$$

Let us note that the three quantities between brackets are positive or null.

**Illustration: back to the normal/standard model of Chaboche-Marquis**

We have seen that the potential $\phi^*(\sigma_D, \tilde{\mathbf{X}}, \tilde{p})$ is equal to:

$$\phi^*(\sigma_D, \tilde{\mathbf{X}}, \tilde{p}) = \frac{k}{n+1} \langle z \rangle_+^{n+1} \quad (106)$$

with

$$z = |\sigma_D - c^{1/2}\tilde{\mathbf{X}}| + \frac{1}{2}ac|\tilde{\mathbf{X}}|^2 - \ell(\tilde{p}) \quad (107)$$

It follows:

$$a(\sigma_D, \tilde{\mathbf{X}}) = \ell^{-1}(|\sigma_D - c^{1/2}\tilde{\mathbf{X}}| + \frac{1}{2}ac|\tilde{\mathbf{X}}|^2) \quad (108)$$

which is convex because $\ell$ is an increasing concave function. One also has:

$$\phi^*(\sigma_D, \tilde{\mathbf{X}}, \tilde{p}) = \gamma(\langle a(\sigma_D, \tilde{\mathbf{X}}) - \tilde{p} \rangle_+) \quad (109)$$

with

$$\gamma(r) = \frac{k}{n+1} [\ell(r)]^{n+1} \quad \forall r \in \mathbf{R}^+ \quad (110)$$
## Optimal solution of standard ill-posed problems

Let us consider that the behavior of the material is known and defined by the potential $\phi^*(\boldsymbol{\sigma}, \tilde{\boldsymbol{X}}, \tilde{\boldsymbol{p}})$ which is a convex function such that:

$$ \phi^* \geq 0; \quad \phi^*(0, \mathbf{0}, \mathbf{0}) = 0 \qquad (111) $$

The displacement $\underline{\mathcal{U}}^K$ and the stress $\boldsymbol{\sigma}^S$ defined over $[0, T] \times \Omega$ should be admissible, i.e.:

$$ \underline{\mathcal{U}}^K \in \mathcal{U}_{ad}^{K[0,T]}; \quad \boldsymbol{\sigma}^S \in \mathcal{S}_{ad}^{S[0,T]} \qquad (112) $$

and the problem to solve is a standard ill-posed problem. A solution could be given by the pair $(\underline{\mathcal{U}}^K, \boldsymbol{\sigma}^S)$ which minimizes the global CRE functional $E_{CRE}$ over the space $\mathcal{U}_{ad}^{K[0,T]} \times \mathcal{S}_{ad}^{S[0,T]}$. $E_{CRE}$ being a strongly convex function for standard material models, the solution is unique. Here, the optimal solution is computed.

First we construct a solution defined by a modified potential $\phi_M^*$ and in a second time, we prove that it is optimal. For elasto-(visco-)plastic materials, the modified constitutive relation is also elasto-(visco-)plastic but nonlocal; however, the optimal solution can be computed without difficulty thanks to an incremental numerical technique.

### A remarkable constraint

We will see that for standard ill-posed problems, $\epsilon_p$ is not a priori arbitrary. One has:

$$ \epsilon_p = \epsilon^K - \mathbf{K}^{-1} \dot{\boldsymbol{\sigma}}^S \qquad (113) $$

and the question is to compute $\epsilon^K$ and $\boldsymbol{\sigma}^S$ from $\epsilon_p$. First, let us compute $\boldsymbol{\sigma}^S$; it is defined for $t \in [0, T]$ by:

$$ \boldsymbol{\sigma}^S \in \mathcal{S}_{ad}^S; \quad \int_{\Omega} \mathrm{Tr}[\boldsymbol{\sigma}^S \mathbf{K}^{-1} \boldsymbol{\sigma}^*] d\Omega = - \int_{\Omega} \mathrm{Tr}[(\epsilon_p - \epsilon_e^K) \boldsymbol{\sigma}^*] d\Omega \quad \forall \boldsymbol{\sigma}^* \in \mathcal{S}_{ad,0}^S \qquad (114) $$

where $\epsilon_e^K$ is the elastic solution for the K-data:

$$ \underline{\mathcal{U}}_d \text{ over } \partial_1 \Omega \times ]0, T[; \quad \underline{\mathcal{F}}_d \text{ over } (\partial_2 \Omega - \partial_{12} \Omega) \times ]0, T[; \quad \underline{f}_d = 0 \qquad (115) $$

It follows:

$$ \epsilon_p + \mathbf{K}^{-1} \boldsymbol{\sigma}^S = \epsilon(\underline{\mathcal{U}}) \qquad (116) $$

with $\underline{\mathcal{U}} \in [\mathbf{H}^1(\Omega)]^3$ and $\underline{\mathcal{U}} = \underline{\mathcal{U}}_d$ over $(\partial_1 \Omega - \partial_{12} \Omega)$.

To get $\underline{\mathcal{U}}^K = \underline{\mathcal{U}}$, $\underline{\mathcal{U}}$ should also satisfy:

$$ \underline{\mathcal{U}} = \underline{\mathcal{U}}_d \text{ over } \partial_{12} \Omega \qquad (117) $$

which is equivalent to:

$$ \int_{\Omega} \mathrm{Tr}[\mathbf{K} \epsilon(\underline{\mathcal{U}}) \epsilon(\underline{\mathcal{W}}^*)] d\Omega = \int_{\partial_{12}\Omega} \mathbf{K} \epsilon(\underline{\mathcal{W}}^*) \underline{n} \cdot \underline{\mathcal{U}}_d dS \quad \forall \underline{\mathcal{W}}^* \in \mathcal{W} \qquad (118) $$

with

$$ \mathcal{W} = \{ \underline{\mathcal{U}} \in [\mathbf{H}^1(\Omega)]^3 \text{ s.t. } \mathrm{div}(\mathbf{K} \epsilon(\underline{\mathcal{U}})) = 0 \text{ over } \Omega, \underline{\mathcal{U}}|_{\partial_1 \Omega - \partial_{12} \Omega} = 0, \mathbf{K} \epsilon(\underline{\mathcal{U}}) \underline{n}|_{\partial_2 \Omega - \partial_{12} \Omega} = 0 \} \qquad (119) $$

Using (116), we get from (118) the following property:

**Property 7** For standard ill-posed problems, the inelastic strain satisfies the constraint:

$$ \int_{\Omega} \mathrm{Tr}[\mathbf{K} \epsilon_p \epsilon(\underline{\mathcal{W}}^*)] d\Omega = \int_{\Omega} \mathrm{Tr}[\mathbf{K} \epsilon_e^{KS} \epsilon(\underline{\mathcal{W}}^*)] d\Omega \quad \forall \underline{\mathcal{W}}^* \in \mathcal{W} \qquad (120) $$

where $\epsilon_e^{KS} = \epsilon_e^K - \epsilon_e^S$, $\epsilon_e^K$ and $\epsilon_e^S$ being the elastic solutions for the K-data and the S-data, respectively.

**Remark 4** For a well-posed computational problem, the constraint (120) disappears.
## Modified potential $\Phi_M^*$

Let us introduce the global potentials:

$$ \Phi^*(\sigma, \tilde{X}, \tilde{p}) = \int_{\Omega} \phi^*(\sigma, \tilde{X}, \tilde{p}) d\Omega; \quad \Phi(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{\tilde{p}}) = \int_{\Omega} \phi(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{\tilde{p}}) d\Omega \quad (121) $$

One considers first the case where the reference potential $\Phi^*$ is regular. The $E_{CRE}$-minimization leads for fixed $(\sigma, \tilde{X}, \tilde{p})$ to the minimization of:

$$ \Phi(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{\tilde{p}}) - \int_{\Omega} [\operatorname{Tr}[\sigma \dot{\epsilon}_p] - \tilde{X} \cdot \dot{\tilde{X}} - \tilde{p} \cdot \dot{\tilde{p}}] d\Omega \quad (122) $$

with the constraint:

$$ \int_{\Omega} \operatorname{Tr} \left[ \mathbf{K}(\dot{\epsilon}_p - \dot{\epsilon}_e^{KS}) \epsilon(W^*) \right] d\Omega = 0 \quad \forall W^* \in \mathcal{W} \quad (123) $$

This constraint being equivalent to $\dot{\epsilon}_p \in H_p$, the previous minimization problem can be rewritten as the minimization of:

$$ \Phi(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{\tilde{p}}) + \chi_{H_p}(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{\tilde{p}}) - \int_{\Omega} [\operatorname{Tr}[\sigma \dot{\epsilon}_p] - \tilde{X} \cdot \dot{\tilde{X}} - \tilde{p} \cdot \dot{\tilde{p}}] d\Omega \quad (124) $$

where $\chi_{H_p}$ is the characteristic function of the subspace $H_p$. It follows that $(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{\tilde{p}})$ is the subdifferential of the Legendre-Fenchel transform of $\Phi + \chi_{H_p}$, for which we know $\Phi^*$ and $\chi_{H_p}^*(\sigma, \tilde{X}, \tilde{p})$ equal to:

$$ \left\{ \begin{array}{ll} \operatorname{Tr}[\sigma \dot{\epsilon}_e^{KS}] & \text{if } \sigma = -\mathbf{K} \epsilon(W) \text{ with } W \in \mathcal{W} \\ +\infty & \text{otherwise} \end{array} \right. \quad (125) $$

This is a composite convex function; from [35], we get that the Legendre-Fenchel transform of $\Phi + \chi_{H_p}$ is:

$$ \inf_{\sigma = \sigma' - \mathbf{K} \epsilon(W), W \in \mathcal{W}} \left[ \Phi^*(\sigma', \tilde{X}, \tilde{p}) + \chi_{H_p}^*( -\mathbf{K} \epsilon(W), \tilde{X}, \tilde{p}) \right] \quad (126) $$

Finally, we have:

$$ \inf_{W \in \mathcal{W}} \left[ \Phi^*(\sigma + \mathbf{K} \epsilon(W), \tilde{X}, \tilde{p}) - \int_{\Omega} \operatorname{Tr}[\mathbf{K} \dot{\epsilon}_e^{KS} \epsilon(W)] d\Omega \right] \quad (127) $$

which defines the modified potential $\Phi_M^*$. In the case of elastoplastic materials, one has a slight modification:

$$ \Phi_M^*(\sigma, \tilde{X}, \tilde{p}) = \Phi^*(\sigma, \tilde{X}, \tilde{p}) + \inf_{W \in \mathcal{W}} \left[ \Phi^*(\sigma + \mathbf{K} \epsilon(W), \tilde{X}, \tilde{p}) - \int_{\Omega} \operatorname{Tr}[\mathbf{K} \dot{\epsilon}_e^{KS} \epsilon(W)] d\Omega \right] \quad (128) $$

## Properties of the solution

The considered solution is characterized by the two dual convex potentials $\phi_M^*$ and $\phi_M$, which defines a nonlocal in space constitutive relation. However, this property does not involve any serious difficulty for the computation. The classical incremental method can be used, the only new point being the calculation at each time step of the $W$-problem for which a reduction method is recommended.

To prove the uniqueness of the solution, let us consider two distinct solutions $(\sigma, \tilde{X}, \tilde{p})$ and $(\sigma', \tilde{X}', \tilde{p}')$. One has:

$$ \begin{aligned} \dot{A} &= \int_{\Omega} [\operatorname{Tr}[(\sigma - \sigma')(\dot{\epsilon}_p - \dot{\epsilon}_p')] - (\tilde{X} - \tilde{X}') \cdot (\dot{\tilde{X}} - \dot{\tilde{X}}') - (\tilde{p} - \tilde{p}') \cdot (\dot{\tilde{p}} - \dot{\tilde{p}}')] d\Omega \\ &= \phi_M' + \phi_M^* - \int_{\Omega} [\operatorname{Tr}[\sigma \dot{\epsilon}_p'] - \tilde{X} \cdot \dot{\tilde{X}}' - \tilde{p} \cdot \dot{\tilde{p}}'] d\Omega \\ &\quad + \phi_M + \phi_M' - \int_{\Omega} [\operatorname{Tr}[\sigma' \dot{\epsilon}_p'] - \tilde{X}' \cdot \dot{\tilde{X}}' - \tilde{p}' \cdot \dot{\tilde{p}}'] d\Omega \ge 0 \end{aligned} \quad (129) $$
Moreover, we have for standard ill-posed problems:

$$ \int_{\Omega} \mathrm{Tr}[(\sigma^S - \sigma^{S'})(\dot{\epsilon}^K - \dot{\epsilon}^{K'})]d\Omega = 0 \quad (130) $$

and thus:

$$ \int_{\Omega} \mathrm{Tr}[(\sigma^S - \sigma^{S'})(\dot{\epsilon}_p - \dot{\epsilon}'_p)]d\Omega = - \int_{\Omega} \mathrm{Tr}[(\sigma^S - \sigma^{S'})(\mathbf{K}^{-1}(\dot{\sigma}^S - \dot{\sigma}^{S'})]d\Omega \quad (131) $$

It follows from (129), (131) and the initial conditions at $t = 0$:

$$ A_{|t} = -\frac{1}{2} \int_{\Omega} \left[ \mathrm{Tr}[(\sigma^S - \sigma^{S'})(\mathbf{K}^{-1}(\sigma^S - \sigma^{S'})) + (\tilde{\mathbf{X}} - \tilde{\mathbf{X}}') \cdot (\tilde{\mathbf{X}} - \tilde{\mathbf{X}}') + (\tilde{\mathbf{p}} - \tilde{\mathbf{p}}')^2] \right] d\Omega \geq 0 \quad (132) $$

Consequently, the two solutions can not be distinct which proves the uniqueness.

Moreover, we have the following main property:

**Property 8** The solution associated to the modified potentials $\Phi_M$ and $\Phi_M^*$ is optimal i.e. it minimizes the CRE.

*Proof* Let $(\underline{\mathcal{U}}^{K'}, \sigma^S)$ be the optimal solution i.e., the minimizer of the global CRE over the admissible space $\mathcal{U}_{ad}^{K,[0,T]} \times \mathcal{S}_{ad}^{S,[0,T]}$ taking into account the constraint (120). One has the following result, where $\Phi$ can be replaced by $\Phi_M$:

$$ a = \int_{0}^{T} \eta(t)dt \left[ \Phi_M(\dot{\epsilon}'_p, -\dot{\mathbf{X}}', -\dot{\mathbf{p}}') + \Phi^*(\sigma', \tilde{\mathbf{X}}', \tilde{\mathbf{p}}') - \int_{\Omega} \left[ \mathrm{Tr}[\sigma' \dot{\epsilon}'_p] - \tilde{\mathbf{X}}' \cdot \dot{\mathbf{X}}' - \tilde{\mathbf{p}}' \dot{\mathbf{p}}' \right] d\Omega \right] \quad (133) $$

The solution $(\underline{\mathcal{U}}^{K}, \sigma^S)$ associated to the modified constitutive relation satisfies:

$$ 0 = \int_{0}^{T} \eta(t)dt \left[ \Phi_M(\dot{\epsilon}_p, -\dot{\mathbf{X}}, -\dot{\mathbf{p}}) + \Phi_M^*(\sigma, \tilde{\mathbf{X}}, \tilde{\mathbf{p}}) - \int_{\Omega} \left[ \mathrm{Tr}[\sigma \dot{\epsilon}_p] - \tilde{\mathbf{X}} \cdot \dot{\mathbf{X}} - \tilde{\mathbf{p}} \dot{\mathbf{p}} \right] d\Omega \right] \quad (134) $$

Let us compute:

$$ \bar{A} = \int_{0}^{T} \eta(t)dt \int_{\Omega} \left[ \mathrm{Tr}[(\sigma - \sigma')(\dot{\epsilon}_p - \dot{\epsilon}'_p)] - (\tilde{\mathbf{X}} - \tilde{\mathbf{X}}') \cdot (\dot{\mathbf{X}} - \dot{\mathbf{X}}') - (\tilde{\mathbf{p}} - \tilde{\mathbf{p}}')(\dot{\mathbf{p}} - \dot{\mathbf{p}}') \right] d\Omega \quad (135) $$

One has from (133) and (134):

$$ \begin{aligned} \bar{A} = \int_{0}^{T} \eta(t)dt & \left[ \Phi_M((\dot{\epsilon}'_p, -\dot{\mathbf{X}}', -\dot{\mathbf{p}}') + \Phi_M^*(\sigma, \tilde{\mathbf{X}}, \tilde{\mathbf{p}}) \right. \\ & \left. - \int_{\Omega} \left[ \mathrm{Tr}[\sigma \dot{\epsilon}'_p] - \tilde{\mathbf{X}} \cdot \dot{\mathbf{X}}' - \tilde{\mathbf{p}} \dot{\mathbf{p}}' \right] d\Omega \right] \\ & + \int_{0}^{T} \eta(t)dt \left[ \Phi_M((\dot{\epsilon}_p, -\dot{\mathbf{X}}, -\dot{\mathbf{p}}) + \Phi^*(\sigma', \tilde{\mathbf{X}}, \tilde{\mathbf{p}}') \right. \\ & \left. - \int_{\Omega} \left[ \mathrm{Tr}[\sigma' \dot{\epsilon}_p] - \tilde{\mathbf{X}}' \cdot \dot{\mathbf{X}} - \tilde{\mathbf{p}}' \dot{\mathbf{p}}' \right] d\Omega - a \right] \end{aligned} \quad (136) $$

and consequently $\bar{A} \geq 0$. Moreover, we get:

$$ \int_{\Omega} \mathrm{Tr}[(\sigma - \sigma')(\dot{\epsilon} - \dot{\epsilon}')]d\Omega = 0 \quad (137) $$
and thus:

$$ \int_{\Omega} \mathrm{Tr}[(\sigma - \sigma')(\dot{\epsilon}_p - \dot{\epsilon}'_p)]d\Omega = - \int_{\Omega} \mathrm{Tr}[(\sigma - \sigma')\mathbf{K}^{-1}(\dot{\sigma} - \dot{\sigma}')]d\Omega \quad (138) $$

It follows with $\eta(t) = 1 - t/T$:

$$ \bar{A} = -\frac{1}{2T} \int_{0}^{T} dt \int_{\Omega} [\mathrm{Tr}[(\sigma - \sigma')\mathbf{K}^{-1}(\dot{\sigma} - \dot{\sigma}')] + (\tilde{\mathbf{X}} - \tilde{\mathbf{X}}') \cdot (\tilde{\mathbf{X}} - \tilde{\mathbf{X}}') + (\tilde{p} - \tilde{p}')^2] d\Omega \quad (139) $$

and thus $\bar{A} = 0$. One concludes that:

$$ (\sigma, \tilde{\mathbf{X}}, \tilde{p}) = (\sigma', \tilde{\mathbf{X}}', \tilde{p}') \quad \text{over } [0, T] \times \Omega \quad (140) $$

Let us remind that $(\dot{\epsilon}_p, -\dot{\tilde{\mathbf{X}}}, -\dot{\tilde{p}})$ is the minimizer of the CRE for fixed $(\sigma, \tilde{\mathbf{X}}, \tilde{p})$. It follows:

$$ (\dot{\epsilon}_p, -\dot{\tilde{\mathbf{X}}}, -\dot{\tilde{p}}) = (\dot{\epsilon}'_p, -\dot{\tilde{\mathbf{X}}}', -\dot{\tilde{p}}') = \partial \Phi_M^*(\sigma, \tilde{\mathbf{X}}, \tilde{p}) \quad (141) $$

□

## Optimal data-driven material model for elastoplastic materials

We only consider the case of elastoplastic materials as elasto-viscoplastic materials can be treated in the same way. The only difference is the identification of one more function, the function $\gamma$.

### Experimental data

Experimental data has been introduced in “Experimental data” section. It is, for $i \in \mathbf{N}_{test}$:

$$ \underline{U}_d(i) \quad \text{over } \partial_1 \Omega \times [0, T]; \quad \mathbf{F}_d(i) \quad \text{over } \partial_2 \Omega \times [0, T] \quad (142) $$

The Hooke tensor has been identified beforehand.

We also have seen previously that for stable elastoplastic materials, the mathematical shape of the most general data-driven model compatible with knowledge from physics and materials science is defined by a unique convex function $a$ belonging to $\mathbf{A}$:

$$ a \in \mathbf{A} \equiv \{ a \text{ s.t. } a \text{ is a convex function over } \mathbf{R}^{6+q}, a(0, \mathbf{0}) \le 0 \} \quad (143) $$

### Computation of the optimal data-driven model

The optimal data-driven material model is the minimizer of the global CRE for the family of tested structures. One first has, for $i \in \mathbf{N}_{test}$:

$$ (\underline{\mathcal{U}}^K(i), \sigma^S(i)) = \underset{(\underline{\mathcal{U}}^K(i), \tilde{\sigma}^S(i)) \in \mathcal{U}_{ad}^K(i) \times \mathcal{S}_{ad}^S(i)}{\arg \min} E_{CRE}^2(q, a; \underline{\mathcal{U}}^K(i), \tilde{\sigma}^S(i)) \quad (144) $$

Then, after defining:

$$ \bar{E}_{CRE}^2(q, a) = \frac{1}{N_{test}} \sum_{i \in \mathbf{N}_{test}} E_{CRE}^2(q, a; \underline{\mathcal{U}}^K(i), \sigma^S(i)) \quad (145) $$

the optimal data-driven material model is recovered as:

$$ (\bar{q}, \bar{a}) = \underset{a \in \mathbf{A}, q \in \mathbf{N}^+}{\arg \inf} \bar{E}_{CRE}^2(q, a) \quad (146) $$

The method to solve this nonlinear problem is similar to the one presented for solving nonlinear elastic material models; the only difference is related to the dimension $q$ of the
vector $\tilde{X}$ of hidden internal variables, which belongs to the unknowns of the optimization problem. Consequently, the optimization problem is first solved for a fixed value of $q$. If the final value of the error is too large, it may mean that the value of the dimension $q$ needs to be increased. The optimal data-driven model should thus be recalculated, leading to a reduction of the error. If this reduction is sufficient, one stops the calculation.

**Remark 5** Finding a convex function $a$ defined over $\mathbb{R}^{6+q}$ is a high-dimensional convex regression problem for which additional assumptions on the function $a$ are more than welcome; they are mandatory to reduce the complexity of the regression problem. Isotropy or orthotropy properties can also be added, but probably the best and most common way is to use micro/meso information through the FE$^2$ method. Of course, the CRE as presented here can also be used to identify parametrized models from materials science. Nevertheless, to be considered data-driven, data should carry significant information in these approaches.

## Data-driven modeling of stable elasto-(visco-)plastic materials taking measurement noise into account

One considers stable materials described by standard models in the situation where the measurement noise should be taken into account. Following [71], the main source of noise comes from full-field displacement measurements. It is noteworthy that in the CRE, this noise impacts only the elastic strain rate $\dot{\epsilon}_e^{KS}$ from the constraint:

$$ \int_{\Omega} \mathrm{Tr}[(\dot{\epsilon}_p - \dot{\epsilon}_e^{KS})\mathbf{K}\epsilon(\underline{W}^*)]d\Omega = 0 \quad \forall \underline{W}^* \in \mathcal{W} \quad (147) $$

where $\dot{\epsilon}_e^{KS}$ depends on the prescribed displacement $\underline{U}_d$. To take the measurement noise into account, we propose to satisfy this constraint approximately through a modified Constitutive Relation Error:

$$ \begin{aligned} mE_{CRE}^2 = & \int_0^T \eta(t) \left[ \Phi(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{p}) + \Phi^*(\sigma, \dot{\tilde{X}}, \tilde{p}) \right. \\ & \left. - \int_{\Omega} \left[ \mathrm{Tr}[\sigma\dot{\epsilon}_p] - \dot{\tilde{X}} \cdot \dot{\tilde{X}} - \tilde{p} \cdot \dot{p} \right] d\Omega \right] dt \\ & + \int_0^T \frac{\mu}{2} ||\dot{\epsilon}_p - \dot{\epsilon}_e^{KS}||^2 dt \end{aligned} \quad (148) $$

where the weight $\mu$ is chosen such that $|||\dot{\epsilon}_p - \dot{\epsilon}_e^{KS}|||$ is of the order of magnitude of the measurement noise. Moreover, one has:

$$ \int_{\Omega} \mathrm{Tr}[\mathbf{K}\epsilon(\underline{W})\epsilon(\underline{W}^*)]d\Omega = |||\dot{\epsilon}_p - \dot{\epsilon}_e^{KS}|||^2 \quad (149) $$

with $\underline{W} \in \mathcal{W}$ and

$$ \int_{\Omega} \mathrm{Tr}[\mathbf{K}\epsilon(\underline{W})\epsilon(\underline{W}^*)]d\Omega = \int_{\Omega} \mathrm{Tr}[(\dot{\epsilon}_p - \dot{\epsilon}_e^{KS})\mathbf{K}\epsilon(\underline{W}^*)]d\Omega \quad \forall \underline{W}^* \in \mathcal{W} \quad (150) $$

Here, the modified CRE (mCRE) replaces the classical CRE developed in the previous sections.

## The hybrid optimal solution

Following the “Modified potential $\Phi_M^*$” section, one computes the dual function of $\Phi_M$ defined by:

$$ \Phi_M(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{p}) = \Phi(\dot{\epsilon}_p, -\dot{\tilde{X}}, -\dot{p}) + \frac{\mu}{2} |||\dot{\epsilon}_p - \dot{\epsilon}_e^{KS}|||^2 \quad (151) $$

First, let us determine the dual function of $\frac{\mu}{2} |||\dot{\epsilon}_p - \dot{\epsilon}_e^{KS}|||^2$; one has:
• if $\sigma = -K\epsilon(W)$ with $W \in W$:

$$g^*(\sigma, \tilde{X}, \tilde{p}) = \int_{\Omega} \left[ \frac{\mu}{2} \mathrm{Tr}[\sigma \mathbf{K}^{-1} \sigma] + \mathrm{Tr}[\sigma \dot{\epsilon}_e^{KS}] \right] d\Omega \quad (152)$$

• otherwise:

$$g^*(\sigma, \tilde{X}, \tilde{p}) = +\infty \quad (153)$$

From [35], we get that the dual function $\Phi_M^*$ is:

$$\inf_{\sigma=\sigma'-K\epsilon(W)} \left[ \Phi^*(\sigma', \tilde{X}, \tilde{p}) + \int_{\Omega} \left[ \frac{1}{2\mu} \mathrm{Tr}[K\epsilon(W)\epsilon(W)] - \mathrm{Tr}[K\epsilon(W)\dot{\epsilon}_e^{KS}] \right] d\Omega \right] \quad (154)$$

One can note that if $\mu = +\infty$, we get the potential associated to the constraint (147). Following the proof of Property 8, we can also prove that the solution which satisfies the modified constitutive relation (potentials $\Phi_M$ and $\Phi_M^*$) is optimal i.e. it is the minimizer of $mE_{CRE}$ over the admissible space $\mathcal{U}_{ad}^{K[0,T]} \times \mathcal{S}_{ad}^{S[0,T]}$.

**Remark 6** The hybrid optimal solution can be useful in nonlinear computational mechanics to compute hybrid numerical twins.

## The optimal data-driven model

All which is described in the “Optimal data-driven material model for elastoplastic materials” section is still valid.

## Conclusion

The CRE has been developed as a general tool to compute data-driven models and also classical parametrized models of materials science seen as minimizers of the global CRE. For stable materials described by a standard model, the normal formulation with only one convex function to identify was put forward. For a given reference constitutive relation, an optimal solution i.e. the admissible displacement-stress pair which minimizes the global CRE has been characterized by its modified constitutive relation, nonlocal in the space variable but easy to compute. Such a solution extended to noisy data leads to the easy computation of a hybrid numerical twin involving elasto-(visco-)plastic materials.

Let us note that materials defined by their microstructure for which a data-driven macro-model should be built is a potential range of application; a great advantage is the large volume of available data that facilitates the convex regression. Based on the homogenized theory of periodic media introduced in [73], this question has already been the subject of numerous works, in particular [20, 56, 62, 77]. Among the extensions of the present work that require further research, we can mention: (i) models of materials defined by a bipotential (not necessarily stable); (ii) large displacement modeling.

## Author contributions

P. L.: Conceptualization, Methodology, Formal analysis, Discussions, Writing L. C.: Conceptualization, Methodology, Formal analysis, Discussions, Writing

## Funding

L. Chamoin acknowledges the funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 101002857).

## Availability of data and materials

Not applicable.
# Declarations

## Competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Received: 30 April 2024 Accepted: 19 November 2024

Published online: 18 December 2024

# References

1. Amos B, Xu L, Kolter J. Input convex neural networks. In: Proceedings of the 34th international conference on machine learning, vol 70. 2017. pp. 146–155.

2. Aquino W, Bonnet M. Analysis of the error in constitutive equation approach for time-harmonic elasticity imaging. SIAM J Appl Math. 2019;79:822–49.

3. As'ad F, Avery P, Farhat C. A mechanics-informed artificial neural network approach in data-driven constitutive modeling. Int J Numer Eng. 2022;123:2738–59.

4. Aubin-Frankowski PC, Szabó Z. Handling hard affine SDP shape constraints in RKHSs. J Mach Learn Res. 2022;23:1–54.

5. Aubin-Frankowski PC, Rudi A. Approximation of optimization problems with constraints through kernel sum-of-squares. Optimization. 2024;1–26.

6. Banerjee B, Walsh TF, Aquino W, Bonnet M. Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional. Comput Methods Appl Mech Eng. 2013;253:60–72.

7. Ben Azzouna M, Feissel P, Villon P. Robust identification of elastic properties using the modified constitutive relation error. Comput Methods Appl Mech Eng. 2015;295:196–218.

8. Benady A, Baranger E, Chamoin L. NN-mCRE: a modified Constitutive Relation Error framework for unsupervised learning of nonlinear state laws with physics-augmented Neural Networks. Int J Numer Methods Eng. 2024;125(8): e7439.

9. Benady A, Baranger E, Chamoin L. Unsupervised learning of history-dependent constitutive material laws with thermodynamically-consistent neural networks in the modified Constitutive Relation Error framework. Comput Methods Appl Mech Eng. 2024;425: 116967.

10. Bonnet M, Constantinescu A. Inverse problems in elasticity. Inverse Probl. 2005;21:R1–50.

11. Chamoin L, Allier PE, Marchand B. Synergies between the Constitutive Relation Error concept and PGD model reduction for simplified V&V procedures. Adv Model Simul Eng Sci. 2016;3:18.

12. Chamoin L, Ladevèze P. Model verification, updating, and selection from the Constitutive Relation Error concept. Adv Appl Mech. 2024;59:311–62.

13. Chouaki AT, Ladevèze P, Prossier L. Updating structural dynamic models with emphasis on the damping properties. AIAA J. 1998;36(6):1094–9.

14. Cueto E, Chinesta F. Thermodynamics of learning physical phenomena. Arch Comput Methods Eng. 2023;30(8):4653–66.

15. Dandin H, Leygue A, Stainier L. Graph-based representation of history-dependent material response in the Data-Driven Computational Mechanics framework. Comput Methods Appl Mech Eng. 2024;419: 116694.

16. Deraemaeker A, Ladevèze P, Romeuf T. Model validation in the presence of uncertain experimental data. Eng Comput. 2004;21:808–33.

17. Diaz M, Charbonnel PE, Chamoin L. Robust energy-based model updating framework for random processes in dynamics: application to shaking-table experiments. Comput Struct. 2022;264: 106746.

18. Diaz M, Charbonnel PE, Chamoin L. A new Kalman filter approach for structural parameter tracking: application to the monitoring of damaging structures tested on shaking tables. Mech Syst Signal Process. 2023;182: 109529.

19. Feissel P, Allix O. Modified constitutive relation error identification strategy for transient dynamics with corrupted data: the elastic case. Comput Methods Appl Mech Eng. 2007;196:1968–83.

20. Feyel F, Chaboche JL. FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. Comput Methods Appl Mech Eng. 2000;183(3–4):309–30.

21. Flaschel M, Kumar S, De Lorenzis L. Automated discovery of generalized standard material models with EUCLID. Comput Methods Appl Mech Eng. 2023;405: 115867.

22. Florentin E, Lubineau G. Identification of the parameters of an elastic material model using the constitutive equation gap method. Comput Mech. 2010;46:521–31.

23. Fuhg J, Bouklas N, Jones R. Learning hyperelastic anisotropy from data via a tensor basis neural network. J Mech Phys Solids. 2022;168: 105022.

24. Gerbaud PW, Néron D, Ladevèze P. Data-driven elasto-(visco)-plasticity involving hidden state variables. Comput Methods Appl Mech Eng. 2022;402: 115394.

25. Geymonat G, Hild F, Pagano S. Identification of elastic parameters by displacement field measurement. Comptes Rendus de l'Académie des Sciences, Mécanique. 2002;330(6):403–8.

26. González D, Chinesta F, Cueto E. Thermodynamically consistent data-driven computational mechanics. Contin Mech Thermodyn. 2019;31:239–53.

27. Grmela M, Öttinger HC. Dynamics and thermodynamics of complex fluids. i. Development of a general formalism. Phys Rev E. 1997;56:6620.

28. Haghighat E, Raissi M, Moure A, Gomez H, Juanes R. A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics. Comput Methods Appl Mech Eng. 2021;379: 113741.

29. Halphen B, Nguyen QS. On the generalized standard materials (in French: Sur les matériaux standards généralisés). Journal de Mécanique. 1975;14(1):39–63.

30. Hannah LA, Dunson DB. Multivariate convex regression with adaptive partitioning. J Mach Learn Res. 2013;14:3261–94.

31. He X, Chen JS. Thermodynamically consistent machine-learned internal state variable approach for data-driven modeling of path-dependent materials. Comput Methods Appl Mech Eng. 2022;402: 115348.
32. Hernández Q, Badías A, González D, Chinesta F, Cueto E. Structure-preserving neural networks. J Comput Phys. 2021;426: 109950.

33. Hild F, Roux S. Digital image correlation. In: Optical methods for solid mechanics. A full-field approach. Wiley; 2012. pp. 183–228.

34. Hildreth C. Point estimates of ordinates of concave functions. J Am Stat Assoc. 1954;49(267):598–619.

35. Hiriart-Urruty JB, Lemarechal C. Convex analysis and minimization algorithms 1. Fundamentals. Berlin: Springer; 1996.

36. Hiriart-Urruty JB. A note on the Legendre–Fenchel transform of convex composite functions. In: Alart P, Maisonneuve O, Rockafellar RT, editors. Nonsmooth mechanics and analysis. Berlin: Springer; 2012. p. 35–46.

37. Huang S, Feissel P, Villon P. Modified constitutive relation error: an identification framework dealing with the reliability of information. Comput Methods Appl Mech Eng. 2016;311:1–17.

38. Huang DZ, Xu K, Farhat C, Darve E. Learning constitutive relations from indirect observations using deep neural networks. J Comput Phys. 2020;416: 109491.

39. Ibanez R, Borzacchiello D, Aguado JV, Abisset-Chavanne E, Cueto E, Ladevèze P, Chinesta F. Data-driven non-linear elasticity: constitutive manifold construction and problem discretization. Comput Mech. 2017;60(5):813–26.

40. Ibanez R, Abisset-Chavanne E, Gonzalez D, Duval J, Cueto E, Chinesta F. Hybrid constitutive modeling: data-driven learning of corrections to plasticity models. Int J Mater Form. 2019;12:717–25.

41. Kirchdoerfer T, Ortiz M. Data-driven computational mechanics. Comput Methods Appl Mech Eng. 2016;304:81–101.

42. Kirchdoerfer T, Ortiz M. Data-driven computing with noisy material data sets. Comput Methods Appl Mech Eng. 2017;326:622–41.

43. Kohn RV, Lowe BD. A variational method for parameter identification. Math Model Numer Anal. 1988;22(1):119–58.

44. Kovachki N, Liu B, Sun X, Zhou H, Bhattacharya K, Ortiz M, Stuart A. Multiscale modeling of materials: computing, data science, uncertainty and goal-oriented optimization. Mech Mater. 2022;165: 104156.

45. Ladevèze P, Leguillon D. Error estimate procedure in the finite element method and application. SIAM J Numer Anal. 1983;20(3):485–509.

46. Ladevèze P, Nedjar D, Reynier M. Updating of finite element models using vibration tests. AIAA J. 1994;32(7):1485–91.

47. Ladevèze P. Nonlinear computational structural mechanics. Berlin: Springer; 1999.

48. Ladevèze P, Pelle JP. Mastering calculations in linear and nonlinear mechanics. New York: Springer; 2005.

49. Ladevèze P, Chamoin L. The Constitutive Relation Error Method: a general verification tool. In: Chamoin L, Diez P, editors. Verifying calculations, forty years on: an overview of classical verification techniques for FEM simulations. Berlin: SpringerBriefs; 2016.

50. Ladevèze P, Néron D, Gerbaud PW. Data-driven computation for history-dependent materials. Comptes Rendus Mécanique. 2019;347(11):831–44.

51. Lambert D, Crouzeix JP, Nguyen VH, Strodiot JJ. Finite convex integration. J Convex Anal. 2004;11(1):131–46.

52. Lemaitre J, Chaboche JL. Mechanics of solid materials. Cambridge: Cambridge University Press; 1990.

53. Leygue A, Coret M, Réthoré J, Stainier L, Verron E. Data-based derivation of material response. Comput Methods Appl Mech Eng. 2018;331:184–96.

54. Linden L, Klein D, Kalina KA, Brummund J, Weeger O, Kästner M. Neural networks meet hyperelasticity: a guide to enforcing physics. J Mech Phys Solids. 2023;179: 105363.

55. Linka K, Kuhl E. A new family of constitutive artificial neural networks towards automated model discovery. Comput Methods Appl Mech Eng. 2023;403: 115731.

56. Liu Z, Wu CT, Koishi M. A deep material network for multiscale topology learning and accelerated nonlinear modeling of heterogeneous materials. Comput Methods Appl Mech Eng. 2019;345:1138–68.

57. Lopez E, Gonzalez D, Aguado JV, Abisset-Chavanne E, Cueto E, Binetruy C, Chinesta F. A manifold learning approach for integrated computational materials engineering. Arch Comput Methods Eng. 2018;25:59–68.

58. Marchand B, Chamoin L, Rey C. Real-time updating of structural mechanics models using Kalman filtering, modified Constitutive Relation Error and Proper Generalized Decomposition. Int J Numer Methods Eng. 2016;107(9):786–810.

59. Marchand B, Chamoin L, Rey C. Parameter identification and model updating in the context of nonlinear mechanical behaviors using a unified formulation of the modified constitutive relation error concept. Comput Methods Appl Mech Eng. 2019;345:1094–113.

60. Masi F, Stefanou I. Multiscale modeling of inelastic materials with thermodynamics-based artificial neural networks (tann). Comput Methods Appl Mech Eng. 2022;398: 115190.

61. Mazumder R, Choudhury A, Iyengar G, Sen B. A computational framework for multivariate convex regression and its variants. J Am Stat Assoc. 2019;114(525):318–31.

62. Michel JC, Suquet P. A model-reduction approach to the micromechanical analysis of polycrystalline materials. Comput Mech. 2016;57:483–508.

63. Mielke A. Formulation of thermoelastic dissipative material behavior using generic. Contin Mech Thermodyn. 2011;23:233–56.

64. Moreau JJ. Convexity and duality. In: Functional analysis and optimization. Academic Press; 1966.

65. Muzellec B, Bach F, Rudi A. Learning PSD-valued functions using kernel sums-of-squares. arXiv:2111.11306.

66. Nguyen H, Chamoin L, Ha Minh C. mCRE-based parameter identification from full-field measurements: consistent framework, integrated version, and extension to nonlinear material behaviors. Comput Methods Appl Mech Eng. 2022;400: 115461.

67. Öttinger HC, Grmela M. Dynamics and thermodynamics of complex fluids. ii. Illustrations of a general formalism. Phys Rev E. 1997;56:6633.

68. Raissi M, Perdikaris P, Karniadakis G. Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys. 2019;378:686–707.

69. Rockafellar RT. Extension of Fenchel's duality theorem for convex functions. Duke Math J. 1966;33(1):81–99.

70. Romero I. Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems. Int J Numer Methods Eng. 2009;79:706–32.

71. Roux S, Hild F. Comprehensive full-field measurements via digital image correlation. In: Silberschmidt V, editor. Comprehensive mechanics of materials. Amsterdam: Elsevier; 2024.
72. Samir Z, Chamoin L, Abbas M. A domain decomposition strategy for mCRE-based model updating in dynamics. Comput Methods Appl Mech Eng. 2023;416: 116348.

73. Sanchez-Palencia E. Non-homogeneous media and vibration theory. New York: Springer; 1980.

74. Simonetto A. Smooth strongly convex regression. 2020. arXiv:2003.00771.

75. Thakolkaran P, Joshi A, Zheng Y, Flaschel M, De Lorenzis L, Kumar S. NN-EUCLID: deep-learning hyperelasticity without stress data. J Mech Phys Solids. 2022;169: 105076.

76. Vlassis NN, Sun W. Sobolev training of thermodynamic-informed neural networks for interpretable elasto-plasticity models with level set hardening. Comput Methods Appl Mech Eng. 2021;377: 113695.

77. Yvonnet J. Computational homogenization of heterogeneous materials with finite elements. Berlin: Springer; 2019.

78. Zhang A, Mohr D. Using neural networks to represent Von Mises plasticity with isotropic hardening. Int J Plast. 2020;132: 102732.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.