domain: spi phys

text:
Jan Niklas Fuhg

Amélie Fau

Nikolaos Bouklas

Michele Marino

Elasto-plasticity with convex model-data

Keywords: Data-driven elasto-plasticity, Physics-informed machine learning, Machine learning, Convex yield function, Anisotropic yield function

research documents, whether they are published or not. The documents may come

Introduction

Machine learning and data-driven techniques are currently one of the main drivers of advancement in engineering science. They allow to make efficient use of available data while allowing to automate training processes aiding quick prototyping. Hence, machine learning solutions have been employed in various fields to build surrogates or reducedorder models for complex physical mappings [START_REF] Amir Barati Farimani | Deep learning the physics of transport phenomena[END_REF][START_REF] Li | Heavy rain image restoration: Integrating physics model and conditional adversarial learning[END_REF][START_REF] Kadeethum | A framework for data-driven solution and parameter estimation of pdes using conditional generative adversarial networks[END_REF]. On the other hand, in computational physics, they are also extensively used as direct solvers for deterministic [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Fuhg | The mixed deep energy method for resolving concentration features in finite strain hyperelasticity[END_REF] or uncertain [START_REF] Zhang | Quantifying total uncertainty in physicsinformed neural networks for solving forward and inverse stochastic problems[END_REF][START_REF] Niklas Fuhg | Interval and fuzzy physics-informed neural networks for uncertain fields[END_REF] partial differential equations.

Recently, machine learning approaches have also seen an increased interest in material modeling. Data-driven constitutive modeling has been used to create surrogates of nonlinear path-independent material laws such as hyperelasticity [START_REF] Fernández | Anisotropic hyperelastic constitutive models for finite deformations combining material theory and data-driven approaches with application to cubic lattice metamaterials[END_REF][START_REF] Niklas | On physics-informed data-driven isotropic and anisotropic constitutive models through probabilistic machine learning and space-filling sampling[END_REF][START_REF] Fuhg | Local approximate gaussian process regression for datadriven constitutive models: development and comparison with neural networks[END_REF]. Alternatively, hybrid models have been proposed where known traditional phenomenological models are locally improved with machine-learning-based correction terms when data is available [START_REF] González | Learning corrections for hyperelastic models from data[END_REF][START_REF] Niklas Fuhg | Modeldata-driven constitutive responses: application to a multiscale computational framework[END_REF][START_REF] Frankel | Machine learning constitutive models of elastomeric foams[END_REF]. In Ref. [START_REF] Niklas Fuhg | Modeldata-driven constitutive responses: application to a multiscale computational framework[END_REF], the term model-data-driven material modeling for these types of hybrid constitutive laws has been established. In contrast to path-independent material behavior where one-to-one mappings between the strain and stress significantly simplify material modeling, path-dependent closure models rely on time-and/or history-evolution of the material. In terms of plasticity, different theoretical frameworks for phenomenological approaches, such as hypo-plasticity [START_REF] Wu | Hypoplasticity then and now[END_REF] or elasto-plasticity [START_REF] Hill | The mathematical theory of plasticity[END_REF], have been developed over the years. Elasto-plasticity, for example, distinguishes between elastic and inelastic constitutive responses by employing yield surfaces. On the other hand, hypo-plastic models make no distinction between elastic and plastic strains and do not use yield functions to characterize the onset of yielding. The divide between these two classical theories in phenomenological elasto-plastic modeling can also be found in their data-driven counterparts.

The earliest machine-learning-based approaches for plasticity (see, e.g., [START_REF] Ghaboussi | Knowledge-based modeling of material behavior with neural networks[END_REF][START_REF] Ghaboussi | New nested adaptive neural networks (nann) for constitutive modeling[END_REF]) treat the material law as a problem defined by sequential stress-strain data, which is in line with the ideas from hypo-plasticity. This idea has also been applied to a majority of the current works. For example, [START_REF] Mozaffar | Deep learning predicts pathdependent plasticity[END_REF] use recurrent neural networks to successfully map between history-dependent stress-strain data points where the internal feedback connections are used to find embedded internal variables implicitly. Other approaches, such as [START_REF] Huang | A machine learning based plasticity model using proper orthogonal decomposition[END_REF] utilize neural networks combined with an experimentally obtainable internal history variable to generate a data-driven plasticity model in two-and three dimensions. However, these approaches suffer from problems with interpretability since it is unclear how the model behaves on unseen data, i.e. data outside the training range. Secondly, efficient sampling schemes are needed that generate representative data that allows to extrapolate to other load cases. This last point has not yet been addressed in a sufficient manner.

Recently, [START_REF] Nikolaos | Sobolev training of thermodynamic-informed neural networks for interpretable elasto-plasticity models with level set hardening[END_REF] propose an approach that is more in line with elasto-plastic theory. Here, the material law is constructed out of modular components, similar to what one would encounter in an analytical formulation, where the elastic and plastic contributions (i.e., yield function) can be trained with separate but complementing predictive tools. This allows for better interpretability and reliability than traditional black-box surrogates. With this idea in mind, [START_REF] Fuhg | Machinelearning convex and texture-dependent macroscopic yield from crystal plasticity simulations[END_REF] introduce a framework to generate texture-dependent yield functions with neural networks which are strictly convex but can also embed microstructural information. Other machine learning-based methods for yield functions have been proposed. A correction term to a yield function was applied in [START_REF] Ibáñez | Hybrid constitutive modeling: data-driven learning of corrections to plasticity models[END_REF]. [START_REF] Hartmaier | Data-oriented constitutive modeling of plasticity in metals[END_REF] explores Support Vector Regression as a tool to determine anisotropic yield functions. Symbolic regression is used in [START_REF] Park | Multiscale constitutive model using data-driven yield function[END_REF] and [START_REF] Gf Bomarito | Development of interpretable, data-driven plasticity models with symbolic regression[END_REF] to find interpretable expressions for isotropic yield functions based on data. [START_REF] Flaschel | Unsupervised discovery of interpretable hyperelastic constitutive laws[END_REF] introduce an approach that fits symbolic regression models to synthetic Digital Image Correlation data. However, all of the discussed models are part of a big-data mindset where a significant amount of data from computer simulations or from full-field measurements [START_REF] Neggers | Big Data in Experimental Mechanics and Model Order Reduction: Today's Challenges and Tomorrow's Opportunities[END_REF] is assumed to be available. Thus, these approaches can not be directly applied using only experimental data that are commonly available from uniaxial and biaxial test as these would probe only a small subset of the stress space.

In this work we propose a hybrid model-data-driven approach that can locally improve phenomenological yield functions where experimental data is available. The framework is developed with data availability, and the idea of data-poor modeling in mind. As a proof-of-concept we propose and study the framework on perfect plasticity. This restriction is chosen deliberately since experimental data describing the shape of the yield function is (typically) only available for the initial yield surface. In future works, we aim to extend it to include hardening.

Typically, it can be seen as a physical requirement that the yield function is convex with regard to its stress-dependent arguments. Hence, in this work, special attention is given to ensure the convexity of the resulting model-data-driven yield function. We recognize that the convexity of the correction term guarantees convexity of the final yield function (assuming that the model component in the hybrid formulation is convex). Hence, the convex extensions to three commonly applied machine learning methods are introduced. These include Input Convex Support Vector Regression, Input Convex Gaussian Process Regression, and Input Convex Neural Networks. The framework is tested for modeling rolled metal sheets which are typically characterized by anisotropic yield. The organization of the paper is as follows. The concept of elasto-plasticity with model-data-driven yield functions in the context of the studied use-case is thoroughly explained in Section 2. The considered convexity constrained machine learning techniques are reviewed in Section 3. The numerical implementation of the framework is detailed in Section 4. Results of model-data-driven elasto-plasticity utilized on rolled sheet metals and compared to pure-data-driven approaches are highlighted in Section 5. The consequences of these findings are discussed in Section 6. The paper is concluded in Section 7.

Materials and Methods

This paper deals with the characterization and modeling of yield functions for metallic materials. The focus is in particular on rolled metal sheets, where x, y and z denote the rolling, transverse, and normal directions (RD, TD, and ND), respectively. These types of materials are generally subjected to a plane stress condition in the (x, y) plane. With obvious notation, σ x , σ y and σ xy respectively represent normal components along the RD and TD, as well as the shear component in the RD-TD plane, of the Cauchy stress tensor σ.

The total strain ε is additively split into an elastic ε e and plastic ε p part:

ε = ε e + ε p . (1)
and the the Cauchy stress is assumed to be linearly dependent on the elastic strain

σ = C : ε e = C : (ε -ε p ) (2)
where C is the elasticity tensor.

Model-data-driven approach

In the proposed model-data-driven approach, plastic flow is described by the following yield function:

f (σ) = f mod (σ) + f rem (σ) ≤ 0 , (3)
where f mod represents a given model component and f rem is the remainder between the exact and model-based yield functions.

In the present work, no hardening effects are accounted. These can be included by adding an explicit dependence on ε p . We restrict ourselves to perfect plasticity, since, data describing the shape of the yield surface is typically only available for the initial surface.

An associative evolution law for the plastic strain is considered, that is

εp = λ R mod σ (σ) + R rem σ (σ) if f (σ) = 0 0 if f (σ) < 0 , (4)
where

R mod σ = ∂f mod ∂σ , R rem σ = ∂f rem ∂σ , (5)
and λ is a scalar (positive) plastic multiplier, whose rate λ is required to be positive, i.e. λ ≥ 0, and different from zero only in presence of plastic flow.

Traditionally, yield functions models are established from uniaxial and biaxial experimental data [START_REF] Banabic | Advances in anisotropy of plastic behaviour and formability of sheet metals[END_REF]. These two are briefly recapitulated in the following.

Experimental data

Addressing uniaxial tension along direction t(θ) = x cos θ + y sin θ (at angle θ with respect to the x-axis in the (x, y) plane), the applied stress tensor results in σ uni θ (s) = s(t(θ) ⊗ t(θ)). The corresponding yield limit, that is stress s such that f (σ uni θ (s)) = 0, is denoted by σ u θ . Another loading case study relevant from the experimental point of view is the application of a biaxial stress state in the RD-TD plane with tan β being the stress ratio between TD and RD. In this case, the stress tensor reads σ bi β (s) = s(x ⊗ x) + s tan β(y ⊗ y) 1 , and the corresponding biaxial yield limit, i.e. stress s such that f (σ bi β (s)) = 0, is denoted by σ b β . Superscripts T and C are employed for indicating tensile (i.e., σ u θ , σ b β > 0) or compressive (i.e., σ u θ , σ b β < 0) stress states, respectively. From the kinematic point of view, plastic flow can be characterized by experimental observations on Lankford coefficients r u θ [START_REF] Banabic | Advances in anisotropy of plastic behaviour and formability of sheet metals[END_REF][START_REF] Fang | A new and direct r-value measurement method of sheet metal based on multi-camera dic system[END_REF]. These coefficients are obtained from uniaxial tension tests along t(θ). Defining t ⊥ (θ) = -sin θx + cos θy as the orthogonal direction to t(θ), Lankford coefficients r u θ = -d p ⊥ /d p z represent the ratio between the in-plane orthogonal d p ⊥ = εp : (t ⊥ ⊗ t ⊥ ) and through-thickness d p z = εp : (z ⊗ z) plastic strain increments, reading

r u θ (σ u θ ) = - d p x sin 2 θ -d p xy sin(2θ) + d p y cos 2 θ d p x + d p y σ=σ u θ (t⊗t) , (6)
with

d p x = εp : (x ⊗ x) = λ ∂f ∂σ x , (7)
d p y = εp : (y ⊗ y) = λ ∂f ∂σ y , (8)
d p xy = εp : (x ⊗ y) = λ ∂f ∂σ xy . ( 9
)
1 In the special case β = 90 • , it results in σ bi 90 (s) = s(y ⊗ y) = σ uni 90 (s).

Model component of the yield function

A general expression of f mod , adopted in what follows, is an invariant-based representation of the form [START_REF] Cazacu | Plasticity-damage couplings: from single crystal to polycrystalline materials[END_REF]

f mod (σ) = 3J 0 2 (σ) -σ2 , (10)
where J 0 2 is the second generalized orthotropic invariant, reading

J 0 2 (σ) = a 1 6 (σ x -σ y ) 2 + a 2 6 (σ y -σ z ) 2 + a 3 6 (σ x -σ z ) 2 + + a 4 σ 2 xy + a 5 σ 2 xz + a 6 σ 2 yz , (11)
with a 1 , . . . , a 6 some material constants.

We use this specific model component since according to [START_REF] Cazacu | Plasticity-damage couplings: from single crystal to polycrystalline materials[END_REF] it represents the "most widely used orthotropic yield criterion for describing the directionality in plastic properties (yield stresses, strains) of sheets and plates".

It has been shown that J 0 2 (σ) is pressure insensitive and so it depends only on the deviatoric stress σ dev = σ -pI with p = Tr(σ)/3, and reduces to its isotropic counterpart J 2 (σ) = Tr(σ 2 dev )/2 for a 1 = . . . = a 6 = 1 (see, e.g. [START_REF] Cazacu | Generalization of drucker's yield criterion to orthotropy[END_REF]). Besides, it has been proven that f mod in eq. ( 10) corresponds to the orthotropic generalization of the von Mises Isotropic Criterion [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF]. In the literature, eq. ( 10) is occasionally reformulated as:

f mod (σ) =F (σ y -σ z ) 2 + G(σ z -σ x ) 2 + H(σ x -σ y ) 2 + + 2Lσ 2 yz + 2M σ 2 xz + 2N σ 2 xy -σ2 , (12)
where F , G, H, L, M and N are material constants (univocally related to the constants a 1 , . . . , a 6 ). Under plane stress conditions (σ z = σ xz = σ yz = 0) only 4 parameters (F , G, H and N ) have an impact which can be experimentally characterized by means of 4 material tests, like for instance:

• RD-uniaxial tensile traction, providing σ uT 0 ; • TD-uniaxial tensile traction, providing σ uT 90 ; • uniaxial tensile traction at θ = 45 • , providing σ uT 45 ; • equi-biaxial tensile traction in the RD-TD plane, with σ x = σ y > 0 and σ xy = 0, providing σ bT 45 .

From these experimental values, the yield function in eq. ( 12) can be straightforwardly calibrated via:

2F = (σ/σ uT 90 ) 2 + (σ/σ bT 45 ) 2 -(σ/σ uT 0 ) 2 , ( 13
) 2G = (σ/σ uT 0 ) 2 + (σ/σ bT 45 ) 2 -(σ/σ uT 90 ) 2 , ( 14
) 2H = (σ/σ uT 0 ) 2 + (σ/σ uT 90 ) 2 -(σ/σ bT 45 ) 2 , (15)
2N = 4(σ/σ uT 45 ) 2 -(σ/σ bT 45 ) 2 , (16)
with σ = σ uT 0 . Alternative strategies based for instance on Lankford coefficients r u θ are available, for which the interested reader is referred to, e.g., [START_REF] Cazacu | Plasticity-damage couplings: from single crystal to polycrystalline materials[END_REF].

Limitations of model component f mod

The modeling component of the yield function f mod (σ) in eq. ( 12) can be effectively applied in a wide range of applications, being, in fact, the most widely used orthotropic yield criterion for describing the directionality in plastic properties of sheets and plates [START_REF] Cazacu | Plasticity-damage couplings: from single crystal to polycrystalline materials[END_REF]. However, the orthotropic von Mises criterion is endowed by intrinsic limitations given by its quadratic form and the relatively low number of available parameters (4 in plane stress), [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF]. Therefore, some typical experimental observations in metals cannot be accounted for, like:

• the fitting of additional uniaxial or biaxial yield stresses from the ones used for calibration, see e.g. Eqs.

(13)-( 16); • the simultaneous fitting of uniaxial yield stress values σ u θ and Lankford coefficients r u θ , since the calibration can be based either on the former or the latter [START_REF] Cazacu | Plasticity-damage couplings: from single crystal to polycrystalline materials[END_REF]; • typical tension-compression asymmetries experienced in material strength e.g. σ uT 0 ̸ = σ uC 0 .

To better explain these aspects, an exemplary case study is presented in Fig. 1. Here, parameters F , G, H and N have been calibrated using known data values at σ uT 0 , σ uT 45 , σ uT 90 and σ bT 45 . It can be seen that the calibration data is perfectly represented, but additional experimental tests cannot be incorporated, resulting in significant inaccuracies due to the adopted yield function f mod (σ). The effect of these inaccuracies on the obtained constitutive responses are shown in Fig. 2 in terms of a strain-driven uniaxial traction test along RD. We remark that, in order to obtain a rigorous comparison between the predicted and true yield surfaces in the forthcoming numerical results, the calibration data in this work is generated by means of an analytical orthotropic yield criterion by Cazacu and Barlat [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF], presented in detail in A.

The traditional approach would be to seek a refined formulation with more parameters. However, the fitting capabilities of the generalized expression might anyway be too low due to intrinsic limitations of the chosen mathematical relationship, and the modeler's choices on the adopted relationship highly affect the outcome. Moreover, the number of parameters is finite, and so will be the data that can be reproduced exactly via a proper calibration step. When more data points than parameters are available, then the parameters are generally calibrated by minimizing an error function, the form and weight factors of which introduce another impact due to decisions of the modeler.

Accordingly, in order to ensure a general accurate description of available information, the authors here propose to introduce a surrogate model for the remainder term f rem in eq. ( 3), built by machine learning. Hence, the remainder term represents the data-driven component in the model-data-driven elasto-plastic approach.

Data-driven component

To introduce the data-driven part of our framework, first we discuss the type of data that could potentially be available from experiments. In general, the data is given as discrete values of the remainder term at given stress states. Hence, our data represents the (negative) effective value of the model yield function at the experimentally measured stress positions.

Three datasets that combine uniaxial and biaxial experiments are generated and contemplated in order to investigate the required and relevant amount of data. From biaxial experiments we get

f rem (σ bi θ * (σ b * θ )) = -f mod (σ bi θ * (σ b * θ )) , * = T, C, (17)
while from their uniaxial counterparts we obtain

f rem (σ uni θ * (σ u * θ )) = -f mod (σ uni θ * (σ u * θ )) , * = T, C. (18)
Combining N independent uniaxial and biaxial experiments, the dataset

D = {x i , f i rem } N i=1 with x i = [σ i xx , σ i yy , σ i xy ]
is hence employed to create the surrogate model frem of the remainder term f rem as detailed later in Section 3. Remark 1. (Using Lankford coefficient data) We remark that there is potential to include Lankford coefficients in the training data set. For the same set of uniaxial tests, known values ruT θ and ruC θ of Lankford coefficients can be determined with:

r rem (σ u * θ ) = ru * θ -r mod (σ u * θ ) , * = T, C (19)
providing an implicit relationship on function derivative of ∂f rem /∂σ at some stress states σ uni θT and σ uni θC . In order to test the capabilities of the model-data approach, we consider a simpler data-regime scenario where no experimental knowledge of the Lankford coefficients is required. This decreases the number of necessary experiments. However, available information about the Lankford coefficients could be straightforwardly included.

Convexity constraints and machine learning approaches

In the middle of the last century Drucker [START_REF] Charles | A more fundamental approach to plastic stress-strain relations[END_REF][START_REF] Charles | A definition of stable inelastic material[END_REF] argued that, based on a stability postulate, the yield surface must be convex, i.e.

(σ ⋆ -σ) : ∂f ∂σ ≤ 0 (20)
where σ ⋆ is an arbitrary macroscopic stress on or inside the yield surface. Hence, ideally, any yield function obtained through the model-data-driven approach of eq. ( 3) should fulfill this criterion. Because addition is a convexitypreserving operation, if both components are convex, the model-data-driven yield function will inherit convexity. It should be highlighted that the model component f mod can directly be chosen to be, as is the case here, a convex function. Hence, the data regression model should aim to fit the best possible convex function to the available data such that convexity of the remainder term f rem is guaranteed. We remark that this is a tighter constraint than strictly necessary (it is indeed a sufficient but not necessary condition, cf. B). Different ways to enforce convexity in machine learning techniques are available, although they require some 12) with parameters obtained from eq. ( 13)) and the true yield (obtained from the orthotropic yield criterion by Cazacu and Barlat [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF], see A). Since the true yield function here follows an analytical definition, the full surface is shown, but it has to be pointed out that realistically the yield function would only be known at a few select points if it was obtained from experimental testing. 12) with parameters obtained from eq. ( 13)) and the true response (obtained from the orthotropic yield criterion by Cazacu and Barlat [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF], see A) during a strain-driven uniaxial traction test: (a) stress-strain relationship, (b) plastic strain.

generalizations for facing the problem at hand. Therefore, we decided to focus on the convexity condition, rather than enforcing less standard constraints on the Hessian of the surrogate model which are out of the scope of the present study.

Guaranteeing the convexity of the surrogate frem is equivalent to ensuring that the Hessian matrix of the prediction for any input x is at least positive semi-definite, i.e., the predictor is subject to the following continuous constraint

H( frem )(x) ⪰ 0 (21)
with

H( frem ) ij (x) = ∂ 2 frem (x) ∂x i ∂x j ( 22
)
where the sign " ⪰ " is used to symbolize positive semidefiniteness. Generally, two distinct concepts can be employed to fulfill this requirement in machine learning frameworks:

• Intrinsically convex regression functions. Some machine learning frameworks allow for a direct fulfillment of the convexity constraint by employing specialized designs.

• Reduction of constraints to finite dimensionality. Instead of ensuring the convexity requirements globally, they are only enforced (with soft or hard constraints) on a finite number of points. Different strategies have been used in the literature to choose the set of points which the requirements should be fulfilled on [START_REF] Wang | Multivariate convex support vector regression with semidefinite programming[END_REF]. For instance, constraints could hold only on training samples [START_REF] Pelckmans | Primal-dual monotone kernel regression[END_REF][START_REF] Du | Nonparametric kernel regression with multiple predictors and multiple shape constraints[END_REF], on sample points in a domain of interest [START_REF] Pensoneault | Nonnegativity-enforced gaussian process regression[END_REF], or on a set of points that was iteratively selected [START_REF] Neumann | Reliable integration of continuous constraints into extreme learning machines[END_REF].

The literature on convex constrained regression models using least-squares estimators is vast, see e.g. [START_REF] Kuosmanen | Representation theorem for convex nonparametric least squares[END_REF][START_REF] Seijo | Nonparametric least squares estimation of a multivariate convex regression function[END_REF][START_REF] Lim | Consistency of multidimensional convex regression[END_REF]. However, these works are mostly concerned with piecewise-linear surfaces which do not allow for the necessary differentiability at the training samples that is required for numerical evaluation of elasto-plastic calculations. In the following, convex extensions to three commonly applied smooth regression techniques are briefly discussed, in particular Support Vector Regression (SVR), Gaussian Process Regression (GPR), and Neural Networks. The first two methods are reliant on the specification of a kernel function r : R d × R d → R which is typically user chosen. The respective kernel matrix R(X, Y ) with X ∈ R n×d and Y ∈ R m×d maps R n×d × R m×d → R m×n and is built with r(X, Y ) i,j = r(X i , Y j ).

Input convex Support Vector Regression

In Support Vector Regression, the value of the regression model at an input x can be obtained by

frem (x) = N i=1 α i r(x, x i ) + b (23)
where both α = [α 1 , . . . , α N ] and b are unknown parameters. The parameters that best fit the observations are typically determined following an ϵ-insensitive loss approach [START_REF] Vapnik | The nature of statistical learning theory[END_REF] where no weight is put onto small residuals. This is typically done by specifying a zone ±ϵ ∈ R + where the residual is not accounted for. The trainable parameters can then be found by solving a linear programming problem [START_REF] Smola | A tutorial on support vector regression[END_REF]. In order to enforce convexity of the SVR predictor in a non-intrinsic manner on a finite number of points {x i c } Nc i=1 , Wang and Ni [START_REF] Wang | Multivariate convex support vector regression with semidefinite programming[END_REF] proposed to reformulate these optimization problems by including the requirement that the Hessian matrix should be positive semi-definite. Since the prediction value of eq. ( 23) is linear this requirement can be included such that the optimization problem takes the form of a semi-definite program

arg min α,b,ξ,a 1 N 1 T ξ + C1 T a s.t. diag (Rα + b1 + ξ -f rem , -Rα -b1 + ξ + f rem , ξ -1ϵ, α -a, -α + a, N i=1 α i F i ) ⪰ 0 ( 24
)
where C is used to tune the trade-off between the empirical error minimization and the regularization term maximization, and where the set of trainable parameters now includes the N -dimensional vectors ξ and α and

F j = diag(H(k(x 1 c , x j ))(x 1 c ), . . . H(k(x Nc c , x j ))(x Nc c )). (25)

Input convex Gaussian process regression

Using some prior assumptions [START_REF] Fuhg | Adaptive surrogate models for parametric studies[END_REF], predictions with Gaussian process regression can be made using the analytical formulation

frem (x) = µ + r T ⋆ (x)R -1 (f rem -1 µ), (26)
with 1 an all-ones vector and µ = (1 T R -1 1) -1 1 T R -1 f rem . Furthermore, the elements of r read r ⋆,i (x) = r(x, x i ) for i = 1, . . . , N.

(

) 27
The kernel function is typically also dependent on a d-dimensional hyperparameter vector θ which needs to be determined based on the available data, e.g. by minimizing the concentrated negative marginal log-likelihood

-ln(f rem |X, θ)) = 1 2 N ln σ 2 (θ) + ln(det R(θ)) (28)
where

σ 2 = 1 N (f rem -1 µ) T R -1 (f rem -1 µ) . (29)
In order to constrain the prediction of eq. ( 26) to be strictly positive it was proposed in [START_REF] Pensoneault | Nonnegativity-enforced gaussian process regression[END_REF] to find the hyperparameters by postulating a constrained optimization problem on the negative marginal log-likelihood over a finite number of points {x i c } Nc i=1 . We reformulate this approach to find the hyperparameters such that the Hessian matrix of the mean prediction value is positive semi-definite on a number of constraining sample points

arg min θ -ln(f rem |X, θ)), s.t. H( frem )(x i c ) ⪰ 0, i = 1, . . . , N c (30)
A noise term can be included as a trainable parameter to improve the flexibility of the GRP regression model, see [START_REF] Edward | Gaussian processes in machine learning[END_REF]. To combat that the noise of the model impacts its accuracy, [START_REF] Pensoneault | Nonnegativity-enforced gaussian process regression[END_REF] add an additional restriction to the constrained optimization model that penalizes noise values that would result in large discrepancies of the trained model and the actual ground truth data. This approach is also followed here. Following [START_REF] Pensoneault | Nonnegativity-enforced gaussian process regression[END_REF] the discrepancy parameter value is set to 0.03.

Other methods for constraining GPR exist but they generally do not guarantee that the posterior is analytically known or are only applicable for simple constraints such as output bounds. For more information and an overview over available techniques we refer to [START_REF] Swiler | A survey of constrained gaussian process regression: Approaches and implementation challenges[END_REF].

Input convex neural networks

In contrast to SVR and GPR, constraining neural networks intrinsically to be convex is possible in a simple and efficient manner to build input convex neural networks (ICNN). Following [START_REF] Amos | Input convex neural networks[END_REF], the standard update formula of neural networks can be rewritten to

z i+1 = g i (W z i z i + W x i x + b i ), (31)
where W z 0 = 0, z 0 = 0 and z k = f rem . The set of trainable parameters is given by {W z 1:k-1 , W x 0:k-1 , b 0:k-1 }. Contrary to traditional feed-forward neural networks, this network architecture includes "passthrough" layers, i.e. the input x is directly connected to the hidden and output layers. These layers are needed to increase the expressiveness of the network since in order to guarantee convexity with regard to the inputs all weights {W z i } k-1 i=1 are required to be non-negative and the activation functions g i are chosen to be non-decreasing and convex. Rectified Linear Unit (ReLU) are the common choice of activation function in this context. The optimal values of the trainable parameters are generally determined by solving an optimization problem over a loss function where for regression problems the mean-squared error is a common choice.

Numerical implementation

The state S of an elasto-plastic material at time t is described by the current values of the total strain, as well as by a set of internal plastic variables S pl (t), here chosen to be plastic-strain, kinematic-hardening stress variables, and scalar plastic multiplier that is S(t) = {ε(t), S pl (t)} with S pl (t) = {ε p (t), λ(t)} .

For a given loading path, updating the material state requires the solution of an evolution problem, which is here faced by means of a time-incremental strategy based on the backward first-order Euler scheme. Accordingly, let us assume that all quantities are known at time t k and denoted with subscript k, and unknowns at the next time step t k+1 = t k + ∆t are indicated without subscripts.

A displacement-driven loading path is addressed. Therefore, the variation of the total strain ε = ε k + ∆ε is given. The evolution of material state is hence fully characterized by finding the value of the following vector of unknowns:

h = vec [{ε p , λ}] , (33)
which collects the values of S pl at time t k+1 .

In the proposed model-data-driven strategy, updating the material state during the evolution of the plastic state, i.e., when S pl (t k+1 ) ̸ = S pl (t k ), requires to solve:

1. the evolution equation (4) for ε p , yielding the following algebraic expression:

Q ε (ε, h) = ε p -ε p,k -(λ -λ k ) R mod σ (σ) + R rem σ (σ) = 0 ; (34)
2. the yield condition:

q λ (ε, h) = f mod (σ) + f rem (σ) = 0 . ( 35
)
The dependency of Q ε and q α on ε and h is highlighted, since the knowledge of these latter quantities allows to uniquely characterize σ from Eqs. [START_REF] Li | Heavy rain image restoration: Integrating physics model and conditional adversarial learning[END_REF]. The evolution of material state is then fully characterized by finding h such that Eqs. ( 34) and ( 35) are satisfied, that is:

Find h such that q pl (ε, h) = vec [{Q ε , q α }] = 0 . ( 36
)
Algorithm 1 Pseudocode for the predictor-corrector scheme.

1: Trial value:

h tr = h k ▷ Initialization 2: if q λ (ε, h tr ) < 0 then ▷ Predictor (elastic step) 3:
return with h = h tr 4: else ▷ Corrector (plastic step) 5:

while ∥q pl (ε, h tr ∥ > tol do

6: ∆h = -[A(ε, h tr )] -1 q pl (ε, h tr ) 7:
h tr ← h tr + ∆h 8:

end while with h = h tr 9: end if

The solution of the non-linear problem [START_REF] Wang | Multivariate convex support vector regression with semidefinite programming[END_REF] for a given value ε is obtained by a predictor-corrector algorithm detailed in Algorithm 1. To this aim, the linearization of the algebraic expression q pl (ε, h) with respect to h is required, i.e.

A(ε, h) = ∂q pl ∂h . (37)
The latter can be obtained by combining partial derivatives of Q ε and q λ with respect to ε p and λ. Most of these relationships follow standard definitions in the field of elasto-plastic numerical models, and are hence not detailed here. Only those modified by the use of the model-data-driven approach are specified here:

∂Q ε ∂ε p = A mod εε I + (λ -λ k ) ∂ 2 f ∂σ∂σ : ∂σ ∂ε p + A rem εε (λ -λ k ) ∂ 2 f rem ∂σ∂σ : ∂σ ∂ε p , (38a)
∂Q ε ∂λ = A mod ελ - ∂f ∂σ A rem ελ - ∂f rem ∂σ , (38b)
noting that ∂σ/∂ε p = -C and ∂q/∂α = -H. In conclusion, the tangent operator for the plastic state evolution can be regarded as the sum of a classical one A mod associated with the model component, and a remainder A rem computed from the first and second derivatives of f rem with respect to stress variables, cf., eq. ( 38). These derivatives can either be obtained as analytical derivatives of the machine learning models by hand or using automatic differentiation techniques. Since all the models used in this work were implemented in Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] we employed automatic differentiation. For more information we refer to [START_REF] Paszke | Automatic differentiation in pytorch[END_REF].

Numerical tests

The proposed approach is tested in this section. The ICNN as well as the unconstrained neural network use ReLU activation functions, are implemented in Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] and the network parameters were optimized using the Adam optimizer [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF]. The network consist of 3 hidden layers with 30 neurons each. The input convex Support Vector Regression models and their unconstrained counterparts where trained using the linear and semi-definite solvers offered by the CVXPY framework [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF]. We used a fifth-order polynomial kernel for all the following results when not specified otherwise [START_REF] Shamshirband | Sensor data fusion by support vector regression methodology-a comparative study[END_REF]. Following Wang and Ni [START_REF] Wang | Multivariate convex support vector regression with semidefinite programming[END_REF] we choose C = 10 -4 and ϵ = 10 -3 .

The constrained and unconstrained Gaussian Process Regression models employ a Matérn 3/2 kernel [START_REF] Fuhg | Adaptive surrogate models for parametric studies[END_REF]. The necessary optimization problems are solved using an interior point method [START_REF] Arkadi | Interior-point methods for optimization[END_REF].

All these hyperparameter and parameter choices were not the results of any elaborate study. Differences between the machine learning techniques presented in the next sections might not be accurately representing their best performances as data-driven predictive tools and should not necessarily be judged as such. Thorough studies would be required to optimize the performance of each machine learning technique, and the number of investigations and case studies to be shown would explode. Such study would hamper to reach the main take-home message of the paper, that is a proof of concept on the integration of model-and data-driven predictions. Clearly, it is however true that the following results can be seen as a broad guidance on possible advantages and disadvantages of the presented machine learning methods by using the techniques for the same tasks and on the same datasets.

We then define three correction datasets for the yield function, c.f. eqs. ( 17) and ( 18):

D N = σ uni θ * (σ u * θ ) σ bi θ * (σ b * θ ) , f rem (σ uni θ * ) f rem (σ bi θ * ) ∀θ∈Θ N (39)
consisting of N = 28, 16 and 8 data points respectively by defining

Θ 28 = [0 • , 30 • , 45 • , 60 • , 90 • , 120 • , 135 • , 150 • ] Θ 16 = [0 • , 45 • , 90 • , 120 • , 150 • ] Θ 8 = [0 • , 45 • , 90 • ]. (40)
Hence, the subscripts on D i and Θ i refer to distinct datasets that contain 28, 16 and 8 unique data points. As for the case study in Section 2.3.1, training data is generated here by means of the orthotropic yield criterion by Cazacu and Barlat [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF] (see A), so as to allow a consistent comparison of the final shape of the yield surfaces. However, this data can alternatively be straightforwardly obtained from well-established testing procedures.

In order to both decrease training time and increase training reliability [START_REF] Goodfellow | Deep learning[END_REF], we employ Min-max feature scaling such that each feature lies inside the range [0, 1].

In the following, 1, 000 randomly sampled points in the three-dimensional domain [σ xx , σ yy , σ xy ] ∈ [-800, 800] 3 MPa were used to enforce convexity in a discrete way for the surrogate modeling techniques that rely on a finite number of points where convexity should be enforced, i.e., SVR and GPR. This particular domain was chosen to strictly include all relevant information of the modeling yield function component (c.f. Figure 1(a)).

First, we study the performances of a series of machine learning techniques in a pure data-driven context for predicting the shape of the true yield function given the experimental data points. Afterwards, we compare these results to the ones obtained using the proposed machine learning techniques for the model-data-driven framework.

Fitting with pure data-driven approaches

In this preliminary study, we take the input data of the 28 point dataset (eq. ( 39)) but assume the output is the true yield function we aim to fit, i.e.

D ⋆ 28 = σ uni θ * (σ u * θ ) σ bi θ * (σ b * θ ) , 0.0 0.0 ∀θ∈Θ28 (41)
where Θ 28 is given in eq. ( 40). Among the machine learning techniques reviewed in Section 3 only Support Vector Regression and Neural Networks allow a reasonable fit of a function with a constant output (due to the hyperparameters of Gaussian Process Regression representing correlation lengths). Hence, GPR will be excluded here. Additionally, the data will be also be fitted using generalized yield function fitting via Fourier series representations. This method is summarized in C. In this method the convexity is also soft-constrained by enforcing the requirement on 1, 000 points in the relevant domain mentioned above. Here, we cut off the series after its first 6 components and solve for the trainable parameters using a stochastic gradient descent scheme (ADAM) and use L 0.25 regularization to generate sparsity in the parameter set, similarly to [START_REF] Flaschel | Discovering plasticity models without stress data[END_REF].

Figure 4 plots the resulting fitted yield functions employing 28 data points of D ⋆ 28 using two unconstrained techniques (an unconstrained Neural Network and unconstrained Support Vector Regression) and three constrained approaches (the input convex Neural Networks formulation, the input convex Support Vector Regression, and the input convex Fourier series approach). It is clear that unconstrained machine learning techniques, see Figs. 4(a) and 4(b), do not generate well-defined predictive tools for yield function representations. They are accurate at predicting zero at the required input points of the dataset, however, due to the missing convexity, the scalar output value fluctuates significantly above and below zero for out-of-sample predictions, their usage appears highly compromised. The advantage of enforcing convexity is visible when looking at the remaining approaches, see Figs. 4(c) and 4(e). The yield surface shapes appear to satisfy convexity. However, we can note that none of the techniques manage to fit the true yield function in an acceptable fashion. The Fourier series representation seems to capture the anisotropy of the yield surface. Among them, the input convex Neural Network appears to generate the most accurate fit of the data points. Results obtained with input convex Neural Network will be used as benchmark performances for data-driven approaches.

Fitting with model-data-driven approach

We employ the three datasets D 28 , D 16 and D 8 in eq. ( 39) to build model-data-driven yield functions based on the model component of eq. ( 12). Moreover, we sequentially test the performance of the three suggested machine learning approaches.

First we discuss the performance of input convex Gaussian Process Regression in the model-data-driven framework.

Figure 5 plots the model-response, the response of the true model and the output of the model-data-driven approach with input convex Gaussian Process Regression in uniaxial and biaxial directions for the three datasets. The respective positions (in this projection) of the input training points are highlighted in each plot. First of all, it can be seen that by locally improving the model component the predicted output of the model-data-driven approach is able to capture the true function. Here, a larger dataset (naturally) helps to improve the accuracy of the predicted output. However, even in the case of 8 (experimental) points the difference between the model-data-driven framework and the true function is still relatively small but in regions where no data is available there are clear deviations. The yield surfaces plotted in three dimensions for Gaussian Process Regression (Figure 6) reveal that a similar pattern can be seen for the whole stress space.

Following, the relevant uniaxial and biaxial curves using input convex Neural Networks for the correction term are compared to the true function and the modeling component in Figure 7 for all three datasets. We can highlight again that the accuracy of the modeling component is clearly improved by adding the data-driven remainder part. With a decreasing number of points the accuracy reduces which is especially noticeable in areas that are further away from the experimental data. The corresponding approximations of the yield surfaces over the full stress field using input convex Neural Networks are displayed in Figure 8. Lastly, the same observations can be made when using input convex Support Vector Regression as the data-driven component, see Figures 9 and10. However, we can see that, compared to the other two approaches, the accuracy of the prediction suffers more with a low amount of data.

In order to highlight the accuracy of the model-data-driven approach in terms of elasto-plastic constitutive responses, Figure 11 plots the stress-strain curves for two different loading paths and for the three different input convex machine learning techniques using the D 28 dataset. It can be seen that, by employing the proposed approach, the corrected curves capture more accurately the true response than the phenomenological model alone. By visual inspection the input convex Neural Network seems to perform slightly better than the other two approaches.

Consider a cycle with ϵ xx oscillating from 0 to +0.1 and -0.1, i.e. load paths as employed in Figures 11 (a-c) as reference. The average residual norm error of the applied Newton-Raphson loop over the number of iterations is plotted in Figure 12. We can see that adding a data-driven component to the phenomenological model decreases the convergence speed and final error of the root finding method. Here, input convex Neural Networks offer a performance that most closely resembles the traditional analytical model.

Finally, Figure 13 plots the Lankford coefficients along the rolling directions for the model, the model-data-driven approaches and the ground truth response. This allows us to investigate how well the model-data-driven framework captures the plastic anisotropy. It can be seen that, except for the low data case test of SVR, all model-data-driven yield functions improve the accuracy of this plastic flow characterization compared to the pure model component. Additionally, it is clearly visible that adding more points increases the quality of the trained yield functions.

Discussion

The results of the previous section show that, for the chosen true response and model component assumptions, the model-data-driven component significantly improves the accuracy of the yield function compared to the pure model counterpart. Additionally, using the proposed approach, the shape of the true yield function is far better approximated in comparison with the purely data-driven techniques (cf., Fig. 4 with Figs. 6, 8 and 10). We can see that the modeldata-driven prediction on 8 training points in all presented methods is significantly more accurate than the data-driven approaches fitted on 28 experimental points.

Convexity constrained Gaussian Process Regression and Neural Network seem to outperform Support Vector Regression specifically when less data is available. However, this might change when a different kernel function is applied. But based on the presented results GPR and NNs have to be the preferred choice. One major advantage of Gaussian Process Regression which has not played a role in this paper is that it allows access to posterior variances. This might be helpful in application cases where the next most effective experiments need to be located in order to improve the accuracy of the prediction based on the currently available information by using for example an active learning approach [START_REF] Fuhg | State-of-the-art and comparative review of adaptive sampling methods for kriging[END_REF]. Furthermore, due to the fact that GPR is deeply rooted in statistical analysis theory, the approach enjoys access to mathematical convergence guarantees. However, as pointed out constraining a GPR model to be convex is a complex task and is only approximate. On the other hand, input convex NNs always fulfill convexity and are more expressive than GPR when a lot of data is available. Even though, using GPR for big data problems has also been explored [START_REF] Fuhg | Local approximate gaussian process regression for datadriven constitutive models: development and comparison with neural networks[END_REF]. However, the focus of the present work is experimental data which from standard measurement setups generally only allow to acquire small, restrictive data sets. Therefore, this common criticism of GPR might not be valid in the context of the present work. We believe that the differences between the convergence results of the Newton-Raphson loop residual stem from the different levels of fulfillment of the convexity requirement. Since, input convex Neural Networks are convex by design they show very good convergence properties when using them in gradient-based solvers. On other hand, Gaussian Process Regression and Support Vector machines are only approximately convex on a finite number of points. However, GPR has (including noise) only 4 trainable parameters whereas SVR has at least 3N + 1 (where N is the number of data points). We speculate, that this in addition to the fact that SVR can be phrased as a simple semi-definite program which allows to exactly fulfill convexity on all constrained points allows SVR to find a function representation that is convex in a wider range of the domain than GPR, which leads to the evident difference in convergence of the residual norm error.

The predicted curves of the Lankford coefficient over the rolling direction angle (Figure 13) highlight that adding experimental values of the Lankford coefficients to the training dataset might even further increase the accuracy of the proposed approach. We will investigate this in future studies.

Overall, the framework seems to yield significant improvements compared to traditional approaches especially in the data-poor regime which is commonly encountered when physical experiments are necessary.

Conclusion

We present a hybrid extension to classical phenomenological model-based elasto-plasticity which allows to locally improve the chosen analytical yield function model where experimental data is available. Our framework is developed with data availability, and the idea of data-poor modeling in mind. The application case studied in this work is related to rolled metal sheets, which exhibit a clear anisotropy in the yield surface which in turn makes fitting accurate yield surfaces more complex. In order to ensure convexity of the predicted yield function in the presented approach, we establish a direct connection to the convexity of the data-driven correction component. Input convex extensions of three classical machine learning techniques are reviewed and applied to the problem at hand. It can be seen that the presented approach significantly improves the accuracy of the standalone model component while achieving better performances than purely data-driven models as well. Furthermore, we are able to generate surprisingly accurate yield surfaces in the data-poor regime.

In future works, the framework can be applied to other use cases which might include a dependence of the yield function on hydrostatic pressure and on actual experimental data. Furthermore, we envision using adaptive sampling techniques in order to determine the best possible next experiment [START_REF] Fuhg | State-of-the-art and comparative review of adaptive sampling methods for kriging[END_REF] that allows to improve the representation of the yield function the most. Additionally, we aim to extend this framework by adding hardening.

Overall, we believe that the presented approach shows a lot of potential to tackle problems in elasto-plasticity modeling where data availability is a concern.

A Orthotropic yield criterion by Cazacu and Barlat [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF] Cazacu and Barlat [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF] propose an anisotropic and asymmetric yield criterion of the form

f true = (J 0 2 ) 3/2 -cJ 0 3 -σ3 , (42)
where c ∈ R, J 0 2 is defined in eq. ( 11) and

J 0 3 = 1
where {a i } 6 i=1 (in J 0 2 ) and {b i } 11 i=1 are two sets of anisotropic coefficients that can be determined when data of experimental yield loci is available. For thin-sheet metals these sets reduce to {a i } 3 i=1 and {b i } 4 i=1 . The authors get them and the additional unknown value c from experimental yield loci data for Mg-4% Li alloys (see Table 1). The yield function estimated by eq. (A) is assumed as the ground truth model in the current study. 1: Anisotropy coefficients and values for c and σ for Mg-4% Li alloy as established in [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF].

B Convexity requirements on yield function

Assume a two-dimensional stress space and define x = [σ x , σ xy , σ y ]. The yield function is convex if its symmetric Hessian matrix H with

H ij = ∂ 2 f ∂x i ∂x j (45)
is at least positive semi-definite which can be denoted by

H(f ) ⪰ 0. ( 46
)
This is equivalent to (see for example [START_REF] Boyd | Convex optimization[END_REF])

a T H(f )a ≥ 0, ∀a ∈ R 3 . ( 47
)
Since the model component is user-chosen, we can pick a convex function f mod . Hence, we know that

a T H(f mod )a ≥ 0, ∀a ∈ R 3 . (48)
In order to ensure the convexity of f , the remainder term has then to fulfill the following condition a T H(f mod )a + a T H(f rem )a ≥ 0 a T H(f rem )a ≥ -a T H(f mod )a [START_REF] Amos | Input convex neural networks[END_REF] for all a ∈ R 3 . To (computationally) simplify this constraint, we tighten the condition on the remainder term to the following instead a T H(f rem )a ≥ 0 ⇐⇒ H(f rem ) ⪰ 0 [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] which includes the condition of eq. ( 49).

Figure 1 :

 1 Figure 1: Comparison of the yield surfaces obtained considering the model f mod (σ) (eq. (12) with parameters obtained from eq. (13)) and the true yield (obtained from the orthotropic yield criterion by Cazacu and Barlat[START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF], see A). Since the true yield function here follows an analytical definition, the full surface is shown, but it has to be pointed out that realistically the yield function would only be known at a few select points if it was obtained from experimental testing.

Figure 2 :

 2 Figure 2: Comparison of the responses obtained along RD only considering the model contribution f mod (σ) of the yield-function (eq. (12) with parameters obtained from eq. (13)) and the true response (obtained from the orthotropic yield criterion by Cazacu and Barlat[START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF], see A) during a strain-driven uniaxial traction test: (a) stress-strain relationship, (b) plastic strain.

Figure 3 :

 3 Figure 3: Visual representation of the network architecture of input convex neural networks.

  (a) Neural Network (b) Support Vector Regression (c) Input convex neural network (d) Input convex Support Vector Regression (e) Input convex Fourier series

Figure 4 :

 4 Figure 4: Yield surface in the three-dimensional stress space using data-driven approaches (a-b) without convexity constraint, (c-d-e) with convexity constraint based on the experimental datasets D ⋆ 28 comprising 28 points, the positions of which are visualized by dots. Comparison with true and model-driven yield functions.

Figure 5 :

 5 Figure 5: Yield function in uniaxial and biaxial directions using the model-data-driven approach with input convex Gaussian Process Regression as the data-driven component based on experimental datasets comprising (a,d) 28, (b,e) 16 and (c,f) 8 points, the positions of which are visualized by dots. Comparison with true and model-driven yield functions.

Figure 6 :

 6 Figure 6: Yield function in the three-dimensional stress space using the model-data-driven approach with input convex Gaussian Process Regression as data-driven component based on experimental datasets comprising (a) 28, (b) 16 and (c) 8 points, the positions of which are visualized by dots. Comparison with true and model-driven yield functions.

Figure 7 :

 7 Figure 7: Yield function in uniaxial and biaxial directions using the model-data-driven approach with input convex Neural Networks as data-driven component based on experimental datasets comprising (a,d) 28, (b,e) 16 and (c,f) 8 points, the positions of which are visualized by dots. Comparison with true and model-driven yield functions.

Figure 8 :Figure 9 :

 89 Figure 8: Yield function in the three-dimensional stress space using the model-data-driven approach with input convex Neural Networks as data-driven component based on experimental datasets comprising (a) 28, (b) 16 and (c) 8 points, the positions of which are visualized by dots. Comparison with true and model-driven yield functions.

Figure 10 :Figure 11 :

 1011 Figure 10: Yield function in the three-dimensional stress space using the model-data-driven approach with input convex Support Vector Regression as data-driven component based on experimental datasets comprising (a) 28, (b) 16 and (c) 8 points, the positions of which are visualized by dots. Comparison with true and model-driven yield functions.

Figure 12 :

 12 Figure 12: Average error norm values over increasing Newton-Raphson-loop iterations employing the model-datadriven approach with three alternative input convex machine learning techniques as data-driven component. Comparison with convergence for the true and model-driven yield functions.

Figure 13 :

 13 Figure 13: Lankford coefficient curves for different machine learning techniques employing the model-data-driven framework with D 8 and D 28 datasets. Comparison with the true response and the model component.

27 (b 1 + b 2 1 + b 4 ) -b 2 -b 3 ] σ 3 zz+ 2b 11 σ xy σ xz σ yz + 2 9 (b 1 + b 2 )σ xx σ yy σ zz - 1 9 (b 1 σ yy + b 2 σ zz )σ 2 xx - 1 9 (b 3 σ 1 9 [(b 1 - 2 zz - σ 2 yz 3 [(b 6 + 2 xz 3 [2b 9 σ 2 xy 3 [

 27121423912929319123623923 zz + b 2 σ xx )σ 2 yy -b 2 + b 4 )σ xx + (b 1 + b 3 + b 4 )σ yy ] σ b 7 )σ xx -b 6 σ yy -b 7 σ zz ] -σ yy -b 8 σ zz -(2b 9 -b 8 )σ xx ] -σ 2b 10 σ zz -b 5 σ yy -(2b 10 -b 5 )σ xx ]

  Consider a split of the yield function into the model and the remainder partf (σ) = f mod (σ) + f rem (σ).

Acknowledgments

M. Marino gratefully acknowledge funding from the Italian Ministry of Education, University and Research MIUR (Programma per Giovani Ricercatori -anno 2017 Rita Levi Montalcini) and of Regione Lazio (POR FESR LAZIO 2014; Progetti di Gruppi di Ricerca 2020; project: BIOPMEAT, n. A0375-2020-36756). N. Bouklas and J.N. Fuhg gratefully acknowledge support by the Air Force Office of Scientific Research under award number FA9550-22-1-0075.

C Fitting of generalized yield function with Fourier series

Recently, there have been attempts in the literature to fit yield surface data based on Fourier series representations [START_REF] Cristian | On the modeling of asymmetric yield functions[END_REF][START_REF] Flaschel | Discovering plasticity models without stress data[END_REF]. The regression model is of the following form

with

and

Here σ 1 ≤ σ 2 ≤ σ 3 are the principal stresses in increasing order and {α i , β i } M i=1 with M ≥ 0 are the trainable parameters. Here M is a user-chosen positive integer value that sets the number of used terms in the approximation.