Over the past decade, advancements in computational frameworks and processing power have made deep neural networks increasingly viable for material modeling. However, purely data-driven models can yield non-physical predictions due to the lack of physical constraints. While recent works that incorporate physics into the training process partially address this issue, they oﬀer no guarantees beyond the scope of the training data. To tackle this challenge, we propose a novel approach that embeds the neural network within the material model, inherently fulﬁlling thermodynamic laws. The networkrepresentsonlytheunknownphysicsallowingustointegrateknowledgeaccumulatedfromdecadesofconstitutive modeling research into the data-driven methodology. By analyzing the trained, embedded networks, we recover existing evolution laws from artiﬁcial training data and discover new evolution laws from experimental data. The discovered evolution laws for isotropic and kinematic hardening can qualitatively predict an experimentally observed yield strength evolution, which conventional evolution laws cannot describe. Keywords: Plasticity, Machine Learning, Neural Network, Data-driven modeling, Equation Discovery 1. Introduction Traditional material modeling requires researchers to propose various model equations and evaluate their abil- itytodescribeobservedmaterialbehavior. Inrecentyears, various alternatives to this approach have been suggested. The model-free approach, introduced by Kirchdoerfer and Ortiz (2016), aims to circumvent the modeling altogether and instead rely directly on an extensive database of ex- perimental data. Another alternative is using a generic model to encode all material behaviors. A neural network is such a generic model according to the universal func- tion approximation theorem (Hornik et al. (1989)). With the rising availability of easy-to-use implementations and computational power, this approach has received increas- ing attention in the last decade. However, the idea is not new, cf. e.g. Ghaboussi and Sidarta (1998), but it also re- quires an extensive database of stress-strain history data. Unfortunately, suchdatabasesarenotreadilyavailabledue to the high cost and time resources needed to run experi- ments. Physics Informed Neural Networks (PINNs) combine physical constraints with available training data by includ- ing known physical laws in the loss functions. For exam- ple, Borkowski et al. (2022) added positive plastic dissi- pation as a regularization term to improve their model. Masi et al. (2021) introduces the “Thermodynamics-based Corresponding author Email address: k.a.meyer@tu-bs.de (Knut Andreas Meyer ANN”, intrinsically fulﬁlling the ﬁrst law of thermodynam- ics and making it easy to regularize wrt. positive dissipa- tion (second law of thermodynamics). Heider et al. (2020) used spectral decomposition to reduce the required train- ing data from experiments and enforce material objectiv- ity (frame invariance), and later used this approach with recurrent neural networks and transfer learning in Fuchs et al. (2021). Malik et al. (2022) also use invariants as inputs to a neural network embedded in their model to describe the yield surface, which was trained on multi- scale simulations of foam structures. Wen et al. (2021) trained a tree-based regressor for plastic slip rate and isotropic hardening modulus before embedding this in a viscoplastic model to evaluate lithium’s temperature and rate-dependent behavior. Although the abovementioned approaches show impres- sive abilities to describe constitutive behaviors, they lead to black-box models. Using a completely diﬀerent idea, Flaschel et al. (2021) introduced a novel method to iden- tify hyperelastic material behavior based on full-ﬁeld data and sparse regression. The concept was recently extended to linear viscoelasticity, Marino et al. (2023). For neural networks, Koeppe et al. (2022) identiﬁed model behavior using a Recurrent Neural Network, and Haghighat et al. (2023) used Feed-Forward PINNs to re-identify material parameters from noisy data. In this paper, we propose to embed a neural network into the standard dissipative constitutive model structure. This idea is similar to Watson and Gupta (1996), who embedded neural networks when modeling microwave cir- ) Preprint April 17, 2023 cuits to represent missing physics. When applied to ma- terial modeling, this approach enables intrinsic fulﬁllment of the laws of thermodynamics and frame invariance, inde- pendent of training data. Furthermore, by describing the evolution equations, we are free to choose any time inte- gration scheme, and in the present study, we use the im- plicit Backward-Euler. Once the model with an embedded neural network has been trained with experimental data, we analyze the neural networks to discover interpretable analytical evolution equations. Finally, we automatically generate new material models that are eﬃcient, accurate, and thermodynamically consistent. Thepaperisorganizedasfollows: InSection2,wederive and motivate the consistent modeling framework that em- beds a neural network. Section 3.1 describes how the ma- terial models are trained to experimental data, while Sec- tion 3.2 outlines the method to discover analytical equa- tions approximating the trained neural network. For the results in Section 4, we ﬁrst show that the proposed ap- proach works for artiﬁcial data generated by a standard J2-plasticity model with Ohno-Wang kinematic hardening in Section 4.1. Thereafter, we apply it to experimental data from Meyer and Ahlström (2023), showing that we can signiﬁcantly improve existing models in Section 4.2. Sections 5 and 6 discuss and summarize our ﬁndings. Fi- nally, in Appendix A, we describe key implementation de- tails for obtaining the parameter sensitivity of material models with embedded neural networks and implicit time integration. 2. Model formulation Strain-driven material models calculate the stress, , for a given time history of strains, . In addition to ﬁt- ting experimental data, it is widely accepted that material models respect the following physical constraints: •A convex Helmholtz’ free energy. •A non-negative dissipation. •Frame invariance. The Helmholtz’s free energy, , for a plasticity formu- lation considering linear elasticity with isotropic and kine- matic hardening, can typically be formulated as (;p;b;k):= e+ p (1) e:=1 2[p] :E: [p] p:=1 3Hkinb:b+1 2jHisojk2 whereandpare the total and plastic strains, respec- tively. Eis the 4th order elastic stiﬀness tensor, while HkinandHisoare the kinematic and isotropic hardening moduli.bandkare internal variables associated with kinematic and isotropic hardening. Thermodynamic sta- bility requires that is a convex function, as discussed inMaugin (1992); Lemaitre and Chaboche (1990); Lubarda (2008). For the speciﬁc form in Equation (1), it is suﬃ- cient that Eis positive deﬁnite and that Hkinis positive to obtain a convex function. The absolute value of Hisoen- ables isotropic softening while maintaining the convexity. This choice is further discussed in Section 2.1. In this work, the isotropic von Mises eﬀective stress is used to deﬁne the yield criterion as :=fvM( )[Y0+] (2) fvM(x):=r 3 2x:xdev whereY0is the initial yield limit. devis the deviatoric part of the tensor , is the back-stress conjugated to b, andis the isotropic hardening stress conjugated to k. The stresses are then deﬁned by the free energy as :=@ @=E: [p] (3a) :=@ @b=2 3Hkinb (3b) :=@ @k=jHisojk (3c) Following these deﬁnitions, the dissipation is given as D=:__ =:_p+ :_b+_k (4) where the Clausius-Duhem dissipation inequality requires thatD 0. Associative evolution laws maximize D, and is typically employed for the plastic strains. For _band _k, this leads to linear hardening laws that do not accurately capture experimentally observed material behaviors, mo- tivating the use of non-associative evolution laws. As an example, the evolution laws for Armstrong-Frederick kine- matic hardening (Frederick and Armstrong (2007)) and Voce isotropic hardening, together with the associative evolution of plastic strains, are given by _p=_ (5a) _b=_gkin() =_ 3 2 1 (5b) _k=_giso() =_ 1 1 (5c) where :=@ @=3 2dev dev fvM( )(6) The functions gkin()andgiso()denote arbitrary evolution laws that are discussed later, while the right hand side shows the speciﬁc Armstrong-Frederick and Voce harden- ing laws. The plastic multiplier, _, is given by the KKT load- ing/unloading conditions, 0;_0;  _= 0 (7) 2 for rate independent behavior. For viscoplastic behavior, _is deﬁned by an overstress function, , _=1 t(; ;) (8) where0is unitless and t>0is a characteristic time. In this work, we adopt the Norton overstress function, =hi Y0n (9) wherenis a material parameter and hi:= max(;0)is the Macaulay bracket. For the evolution laws in Equation (5), the dissipation inequality is fulﬁlled, D=_" Y0+3 2 : 1+2 1# 0 (10) if 1>0and1>0. We will denote this com- plete model, with the Norton overstress function, as the “Chaboche”-model. It is well known that this model can- not accurately describe the behavior of many materials. A critical modeling choice is the form and parameterization of the hardening evolution laws (5b) and (5c). As these functions are unknown, a universal function approxima- tor, such as a neural network, is suitable to model these. However, such a material model must fulﬁll the previously mentioned criteria. Convexity is already assured by not modifying the free energy. For dissipation, the forms in Equations (5) mo- tivates us to propose the following ansatz for evolution equations _b=_h 1[: ] NN k; NNk; i (11a) _k=_[1NNiso] (11b) where NNk;,NNk; , and NNisoare scalar outputs from a neural network. The scalar factors, : and, in front ofNNk;andNNiso, makes it possible to inherently fulﬁll the Clausius-Duhem dissipation inequality, D 0, D=_h Y0+ [: ]2NNk; + [ : ]NN k; +2NNisoi (12) by requiring positive neural network outputs. A non- negative activation function in the last layer fulﬁlls this requirement; see Section 2.3. Furthermore, those scalar factors cause HkinandHisoto retain their physical inter- pretation wrt. the initial plastic stiﬀness. Finally, we require a frame invariant model formulation, implying that the model response is independent of the chosen coordinate system. In particular, tensor compo- nents are neither inputs nor outputs of the neural network; inputs are independent invariants of the current material state, and outputs are the scalar factors in Equation (11). In the most general case for the proposed model, we canthen use the following 19 invariant inputs (see e.g. Boehler (1977))  tr(2 p) tr( 2) tr( ) tr(p) tr( p) tr( 22 p) tr(3) tr(3 p) tr( 3) tr( 2p) tr( 2) tr(2 p) tr(  p) tr( 2 p) tr(2p) tr(2 ) tr(2 2) tr(22 p)(13) noting that tr(2) = 3=2is not used as it is constant. Furthermore, the von Mises yield criterion implies that tr() = tr(p) = tr( ) = 0and that these tensors are symmetric. The choice of using instead of, e.g., dev, is motivated by the dependence between , other invariants, and the eﬀective stress for rate-independent loading. The ﬁnal model uses a reduced set of the above invariants, as discussed in Section 2.3. While the presented model is initially isotropic, struc- ture tensors, cf. e.g. Holzapfel and Ogden (2010), can ex- tend the model to anisotropic responses. For example, a constant anisotropic elastic stiﬀness, E, would add to the invariants by considering, e.g., :E: . However, in the case of isotropy, this reduces to 2Gtr( 2)(whereGis the elastic shear modulus) and does not add an independent invariant. Similarly, an evolving 4th order anisotropy ten- sor,C, can be used in the yield criterion (cf. Meyer and Menzel (2021)). This tensor would be another state vari- able giving rise to additional invariants, such as :C: . 2.1. Isotropic softening The discussion so far has considered an initially harden- ing behavior. However, many metallic materials, such as e.g. carbon steel, may exhibit initial isotropic softening. A negative hardening modulus, Hiso<0, results in isotropic softening, neglecting the absolute value in Equation (1). However, without the absolute value, the Helmholz free energy becomes concave for a negative hardening modu- lus. Therefore, for Hiso<0, we use the following evolution law instead: _k=_ 1 +NNiso+ Y0 (14) In this case, 0, and the evolution of khas the opposite sign in comparison to in Equation (11b). The additional term,=Y 0ensures that the dissipation, D=_"p Y0+pY02 + NN k;[: ]2 + NN k; : +2NNiso# (15) is non-negative for non-negative outputs from the neural network. For a zero-valued network, this corresponds to a isotropic saturation stress of Y0. 3 2.2. Convergence requirements The neural networks are chosen for their ability to ap- proximate any function. However, during training, this implies that the predicted material response can become non-convergent. In practice, we have found it necessary to i) avoid a zero-sized yield surface and ii) avoid too rapid softening. Both these lead to lack of solutions to the lo- calplasticityresidualequations. Theconstraintsdescribed below aid convergence during training but should typically not aﬀect the response of trained models. 2.2.1. Minimum yield strength An isotropic yield limit Y0+ <0yields > 0. This result contradicts the KKT conditions and prevents solu- tions to the local plasticity problem. To ensure a positive yield limit, we constrain the isotropic hardening stress, , via the smooth bounds function hdeﬁned as h(x;a;b ):=(xxb ab1 [ab] exp xa ab +bxb ab<1(16) Fora>b,hprovides a smooth lower bound of xsuch that h(x;a;b ) =x xa b x!1(17) where we note that dh=dx!1asx!a. Speciﬁcally, we constrainas ^=h(+Y0;Ylow;Ymin)Y0 (18) and use ^in the yield criterion such that Y0+ ^Ymin. In all cases, Ylow= 2Ymin= 40 MPa . 2.2.2. Maximum softening In a more detailed analysis compared to Meyer (2020), we analyze the existence of a solution for a rate- independent material. Speciﬁcally, we consider a loading with _= _n, where _>0andnis a 2nd order, normalized tensor, describing the loading direction. The KKT condi- tions in Equation (7), implies that _= 0during plastic loading, resulting in _dev= 2G_ndev2G_ (19) _= 0 = 2G_:ndev _ 3G+2 3Hkin:gkin() +jHisojgiso() (20) where we have isotropic elasticity with shear modulus G, andgkin()andgiso()are deﬁned in Equation (5), with speciﬁc forms in Equations (11) and (14). During plastic loading,:ndev>0, implying that 3G+ [2Hkin=3]:gkin() +jHisojgiso()>0(21)is required as _0. To fulﬁll this requirement, we con- strain the evolution laws by using the smooth bounding functionhfrom Equation (16) as ^giso() =h(jHisojgiso();1:8G;2:0G) jHisoj(22) :^gkin() =h 2Hkin 3:gkin();G 5;G 2 2Hkin=3(23) These bounds ensure that 3G+ [2Hkin=3]:gkin() +jHisojgiso()>0:5G(24) As most materials only exhibit isotropic softening, more isotropic than kinematic softening is allowed. While Equa- tion (22) directly modiﬁes the isotropic hardening law, the limit for the kinematic evolution is implicit. A scaling fac- tor,s, for the neural network outputs yields, ^gkin() =s[: ] NN k;sNNk; (25) s=Hkinh2 3Hkin:gkin();G 5;G 2 Hkin2 3Hkin:gkin()(26) such thats= 1if[2=3]Hkin:gkin()>G=5. Using the modiﬁed functions, ^gisoand^gkinin the evolution of and , ensures that a solution can be found. 2.3. Network architecture and scaling Considered as a function, a feed-forward neural network accepts a vector of inputs, x0, and gives a vector of out- puts,xN, by passing through its Nlayers. The layer iis described as xi=ai:(Aixi1+bi) (27) where the notation a:(x)implies application of the activa- tion function aon each component of the vector xindivid- ually. The trainable parameters of the network are then the weights in the matrices Aiand the bias vectors bifor i2[1;N]. In this study, we reduce the number of invariants from Equation (13) and consider the following six: I1=; I 2= : ; I 3=: ; I4= :p; I 5=:p; I 6=p:p(28) as network inputs. These choices are motivated by the von Mises yield surface. Duringthedevelopmentofthemodelwetestedafewdif- ferent network dimensions and activation functions, and, while many gave good results, a network consisting of ﬁve layers with widths 6-6-6-6-5-3 produced the most stable results. The activation function for the ﬁrst four layers was chosen as ai(x) = tanh(x), while the last activation function,a5(x) =x2ensured the previously discussed non- negativity of the output. 4 In the evolution laws, Equations (11) and (14), the unit of the network outputs is the inverse of stress. To make the networks unit independent and to ensure weights with reasonable magnitudes, the network outputs are divided by a constant, NN0= 1000 MPa . Similarly, many of the inputs also have the unit of stress. These are divided by NNn 0, wherenis an exponent depending on the input’s unit (e.g.n= 2fortr( 2)). 3. Training Section 3.2 Section 3.1 Neural Network model formulation Experimental stress -strain data Chaboche model formulationTrained Neural Network model Trained Chaboche modelInvariant time history data Discovered model formulationTrained discovered modelSparse regressionSection 2 Figure 1: Training methodology workﬂow The overall training procedure in this paper is described in Figure 1, along with the relevant sections. The start- ing points are the Neural Network and Chaboche models, formulated in Section 2. Independent of the material for- mulation, we use the same training process as described in Section 3.1. Section 3.2 outlines the methodology for dis- covering new analytical evolution equations based on the trained Neural Network model. After equation discovery, a new model is automatically generated and trained with experimental data following the process in Section 3.1. 3.1. Training of material models Parameter identiﬁcation for classical material modeling is closely related to training in machine learning: the end goal is to minimize some objective or loss function wrt. a set of parameters. A diﬀerence between regular material models and neural networks is that neural networks are overparameterized such that multiple equivalent minima may exist. Because our models may include overparam- eterized neural networks, we use training algorithms and procedures from machine learning. It turns out that these approaches perform well for regular material models, and we use the same training procedure irrespective of the ma- terial model type. 3.1.1. Parameter bounds and penalization The elastic parameters are ﬁxed to Young’s modulus, E= 210 GPa , and Poisson’s ratio, = 0:3. For some parameters participating in the parameter identiﬁcation, reasonable bounds may be given a-priori. One example is the initial yield strength, Y0. We use gradient-based optimization in the present work, which requires diﬀeren- tiable constraints. Using the smooth bounding functionin Equation (16) for both upper and lower bounds en- sures diﬀerentiability. One issue with this function is that the gradient vanishes when the value of a parameter ap- proaches the hard bounds. To remedy this, we ﬁrst scale the data such that values within the smooth bounds are between 0 and 1, and use the so-called “epsilon insensitive loss”, cf. Lee et al. (2005), to apply L2 regularization on the distance of scaled values to this range (i.e. both -0.1 and 1.1 are penalized by 0:12). Another common issue with neural networks is unstable training due to exploding gradients, which L2 regulariza- tion also remedies Pascanu et al. (2013). In the present study, however, the regular L2 penalization of the neural network’s weights and biases reduced the training perfor- mance. Therefore, all network parameters were treated as regular material parameters, with soft bounds -1 and 1. The only exceptions were the biases in the last layer, where bounds yielding reasonable saturation levels were used (cf. parameters 1and 1in the Chaboche model). In summary, given a design variable x, the actual pa- rameterpused in the model is calculated as gp(p) =ppmin pmaxpmin ^p=8 >>< >>:h x;1;g ph max x>1 x 0x1 h x;0;g ph min x<0(29) p= ^p[pmaxpmin] +pmin wherepminandpmaxare the soft bounds (cf. ain Equa- tion (16)) and ph minandph maxthe hard bounds (cf. bin Equation (16)). Consequently, each parameter pihas a corresponding design variable xi, which is between 0 and 1 if inside the soft bounds. Values outside these bounds are added to the regularization loss, Lr, Lr:=NpX i=1gx(xi)2; gx(x) =8 < :x x< 0 0 0x1 x1x>1(30) 3.1.2. Loss function To have comparable values for the errors between each experiment, simulated stress components, s ij, and experi- mentally measured stress components, e ijare normalized by the number of simulation steps, Nstpand the stress range in the experiment. ^ ij= ijpNstp max e ij min e ij(31) foreither s(simulated) or e(experiment). Using these normalized stresses, the following loss, Lk 5 is calculated for each simulation, k, Lk:=X ij"NX l=1 ^s ij(tl+10)^e ij(tl)2 +N10X l=1hh ^s ij l10 ^s ij li h ^e ij l10 ^e ij lii2# (32) where ^ ij ldenotes the stress in time step l. The ﬁrst term in Equation (32) is the standard L2-loss function. WhileVlassisandSun(2021)showedthatSobolevtrain- ing is advantageous for learning accurate models, the ex- perimentaltrainingdatainthepresentstudydoesnothave the required derivative information. However, ﬁnite diﬀer- ences can approximate the directional derivative along the loading direction, which the second term in Equation (32) uses to regularize the loss function. This addition serves two purposes: Firstly, solutions with oscillation stress re- sponses may occur due to the neural networks’ high ﬂexi- bility. The regularization penalizes such oscillations. Sec- ondly, most experimental scatter in stress-strain curves is variation between test samples. While the curves for re- peatedexperimentshavesimilarshapeswithlownoise, the stress levels may vary slightly. Penalizing the diﬀerence in stress changes gives more importance to the shape of the curve, making the loss function less sensitive to sample variations. For a given batch of Mexperiments, the total loss is then calculated as L:=Lr+1p MMX k=1Lk (33) whereLris deﬁned in Equation (30). 3.1.3. Optimizer and training procedure Three diﬀerent gradient-based optimizers, available in Flux.jl (Innes (2018)), were tested: RMSProp (Tieleman and Hinton (2012)), Adam(Kingma and Ba (2015)), and RAdam(Liu et al. (2020)). For the speciﬁc problem, the RAdamoptimizer gave the most reliable results. While RMSProp and Adamperformed similarly, the ﬁnal losses were somewhat higher for some initial guesses. A coarse sweep revealed that while ﬁne-tuning could improve the convergence rate, the default values worked suﬃciently well for the present study. While eﬃcient implementations of parameter sensitivity for neural networks are abundant, this is not true for nonlinear material models with im- plicit time integration. Appendix A outlines our eﬃcient method for obtaining this sensitivity that the gradient- based optimizers require. The expressiveness of the embedded neural network im- plies that multiple local minima exist. Several strategies for avoiding getting stuck in such local minima exist in the literature, such as “Adaptive activation functions” Jagtap et al. (2020) and mini-batch training Keskar et al. (2017).In this study, we experienced that mini-batch training sig- niﬁcantly improved the performance when we used one experiment in each batch. In addition, one experiment is used as a validation to prevent overﬁtting. The overall training procedure consists of the following steps to learn the scaled parameters x: •ForNminibatch updates: –Simulate and update xfor each mini-batch set. –Simulate the validation set: Save loss. –Simulate the full set: Save loss and update x. •ForNfullbatch updates: –Simulate the validation set: Save loss. –Simulate the full set: Save loss and update x. where the full set typically is the collection of all mini- batches. The experiment in the validation set is never used to update the parameters. 3.2. Equation discovery The equation discovery aims to identify analytical ex- pressions that describe the evolution laws, i.e. the neu- ral network outputs, NNiso,NNk; , and NNk;. Various techniques exist, such as symbolic regression (cf. Holland (1992) (originally published 1975)) and sparse regression (cf. e.g. Brunton et al. (2016)). Herein, we adopted a mod- iﬁed sparse regression technique based on Flaschel et al. (2021), described in more detail below. We start up with a database of candidate functions. While in theory, one could ﬁt the material response to all possible function combinations, this is infeasible due to the high number of possible combinations. Instead, we use our trained neural network model to generate train- ing data for the discovery, consisting of the outputs of the neural network: [NN iso;NNk; ;NNk;]for many diﬀerent combinations of the input invariants: [I1; :::;I 6]. Least square regression can then be used to ﬁnd coeﬃcients for each candidate function. However, this would lead to long expressions with many paramaters, motivating the need for sparse regression techniques. 3.2.1. Generation of input data Assuming that the neural network is fully trained and accurately describes the evolution laws for all possible sce- narios, we could use a large set of random invariant values as input data to generate a database of the network’s out- puts. In practice, however, the network is only trained on invariants encountered during training. Consequently, we run the full training batch while calculating and saving the invariants,I1;:::;I 6in each time step. We then calculate the network outputs corresponding to invariants in each time step. Each invariant, as well as each network output, is scaled by its maximum absolute value. The scaled invariants 6 fulﬁll ^Ii2[1;1]for all time steps. Due to the activa- tion function in the last layer, the scaled outputs fulﬁll ^y2[0;1]in all time steps (where represents iso,k;, ork; ). Finally, duplicate values are removed from the training set. 3.2.2. Generation of candidate functions Each candidate function, fi(^I), should take as input a vector ^Iof scaled invariants (discussed above) and return a scalar output. The approximation of the scaled total output, ^y, is then ^yX i^pifi(^I) (34) For this approach to be suitable, considering the scaling of the invariants, we require that fare positively homo- geneous functions. This implies that polynomials are al- lowed, butnot, e.g., trigonometricfunctionsorexponential functions. Given these restrictions, we build a basis with the following elementary basis functions: •Polynomials of Ijup to 2nd order. •Polynomials of Ij up to 2nd order. •Polynomials of Ij up to 2nd order. •Combinationsof Ij nwithrationalandintegerexpo- nents,n2f0;1;2;1=2;1=3g(max 3 factors per term). •Combinations of Ij nwith rational and integer ex- ponents,n2 f0;1;2;1=2;1=3g(max 3 factors per term). where the use of the Macaulay bracket, hi, enables dif- ferent parameters for positive and negative values of the same inputs. This gives in total 2680 possible functions after ﬁltering out some known equivalent bases (e.g. we know that I3= : 0so all terms containing hI3i can be removed). Additionally, from each group of func- tions with matching values for every sample in the dataset, only the one with the lowest number of active invariants is retained. The limitation of positive homogeneous function can be lifted by not scaling the invariant inputs, which could al- low more general candidate functions, such as exponential functions. In these cases, it often makes sense to include a parameter inside the functions, i.e. exp(p), wherepis a parameter. While this is not directly compatible with the standard regression technique, Marino et al. (2023) use a large set of parameterized exponential functions in combination with k-means clustering to ﬁnd suitable val- ues. Similar to symbolic regression, this approach makes it possible to use more general candidate functions.3.2.3. Discovery Given a number of candidate functions, the goal is to identify which functions can best represent the neural net- work outputs. To study the trade-oﬀ between complexity and accuracy, we would like to construct the Pareto front by comparing the loss and the number of active functions. Foralownumber, NA, ofactivefunctionsandareasonable number,NC, of candidate functions, it is possible to try all possible combinations. The number of combinations is the binomial; NCchooseNA. Consequently, the num- ber of possible combinations increases quickly with both the number of active functions and the available candidate functions. For the experimental data discussed later, af- ter ﬁltering based on the dataset, 869 candidate functions remained. For this case, there are about 110 million pos- sible combinations for three terms, which takes about 10 minutes to evaluate (on a regular laptop) for each network output. Increasing to four terms yields about 24 billion candidates, which, while still feasible, motivates the use of more eﬃcient algorithms. Therefore, we employ the sparse regression algorithm described in Flaschel et al. (2021) to reduce the number of candidate functions. Speciﬁcally, running the sparse regression with multiple initial guesses determines the most frequently selected terms. We then run a binomial search using these terms as candidate func- tions. This selection method is validated by running the binomial search directly for up to three terms. 4. Results 4.1. Artiﬁcial data To validate the proposed workﬂow, we use artiﬁcial data generated by a known plasticity model. Speciﬁcally, the Ohno-Wang kinematic hardening law Ohno and Wang (1993), _b=_ 3 2 fvM( )h: i 1fvM( ) 1m (35) In contrast to the Armstrong-Frederick hardening law, the neuralnetworkcannotexactlyrepresentthisequation. For the training data generation, we use kinematic hardening parameters Hkin= 500 GPa , 1= 500 MPa , andm= 2. In addition, we have elastic parameters E= 210 GPa and= 0:3, initial yield limit Y0= 350 MPa , Voce-type isotropic hardening (Equation (5c)) with Hiso= 25 GPa , and1= 100 MPa , and Norton overstress (Equations (8) and (9)) with t= 1 sandn= 2. 4.1.1. Model training The artiﬁcial training data is generated by simulat- ing strain-controlled biaxial experiments, corresponding to axial-torsion loading of tubular specimens. The load- ing is deﬁned by 11and12, while13=23= 0and 22=33= 0. Four diﬀerent time histories of 11and 12are simulated to create three mini-batch sets for train- ing and one set for validation. We then follow the training 7 0 2500 5000 7500 1000010−610−410−2 EpochLossTraining ValidationFigure 2: Loss function evolution for model training on artiﬁcial data strategyoutlinedinSection3, resultinginthelossfunction evolution in Figure 2. In this case, the lowest validation loss occurs during the mini-batch training phase in epoch 2845. The low relative losses in Figure 2 show that the neu- ral network evolution can replicate the material response for the Ohno-Wang evolution law. Furthermore, Figure 3a shows the axial and shear stress responses compared to the ﬁrst training batch with artiﬁcial experimental data. The results are similar for the two other training batches. As for the training data, no visual diﬀerence can be ob- served between the neural network model and the artiﬁcial experiment data in the validation case in Figure 3b. 0 2 4 6 8 10−1000−50005001000 Time [ s]Stress [ MPa]Exp Sim ( σ11)Sim ( σ12) a) First training case. 0 2 4 6 8 10−1000−50005001000 Time [ s]Stress [ MPa] b) Validation case. Figure3: Simulated(Sim)axial, 11, andshear, 12, stressresponses compared to the artiﬁcial experimental data (Exp).4.1.2. Equation discovery The previous subsection showed that the neural network could approximate the Ohno-Wang evolution law consid- ering the overall material response. The next question is whether this evolution law can be re-identiﬁed from the trainedneuralnetwork, assumingthecorrectequationsare part of the candidate functions. For the speciﬁc material parameters used to generate the artiﬁcial data, in particu- larm= 2, and considering the structure of the kinematic evolution law in Equation (11a), we expect to identify NNiso=p1 NNk; =p2hI3ip I2; I 3=: ; I2= : (36) NNk;= 0 wherep1andp2are constants that depend on the scaling used in the sparse regression. Since we only expect to ﬁnd one term per network output, we just run the binomial search for up to two terms. Table 1 also includes the loss when using zero terms (expected for NNk;). The reported losses are the L2-losses divided by the number of samples and the square of the scaling factor for the network out- put. This scaling enables comparisons with the losses for experimental data in the next section. For the NNisooutput, the loss reduces four orders of magnitude by using one constant term, which is also the expected ﬁnal expression. Adding a second term does re- duce the loss further, but only by a factor of 3.5. The loss also drops four orders of magnitude when adding a single term to the NNk; output. The added term is equivalent to that in Equation (36), noting that I2= : 0, such thathI2i=I2. Finally, the loss for the NNk;output is very low without any terms. Adding one or two terms reduces the loss only marginally. These ﬁndings show that the expressions in Equation (36) have been identiﬁed. Table 1: Discovered equations for artiﬁcial training data loss equation NNiso 1:2102- 2:4102p1 6:9103p1+p2hI1i23p hI2ihI3i2 NNk; 7:1100- 9:0104p1p hI2ihI3i 6:1104p1p hI2ihI3ip2p hI3i NNk; 1:7104- 1:2104p1hI5i 2:1105p13p hI1i3p hI2ip hI3i+p23p hI1i3p hI3i 8 4.2. Experimental data Having demonstrated that the proposed procedure works for artiﬁcial data with a known solution, we apply it to uniaxial experimental data from Meyer and Ahlström (2023). 4.2.1. Model training Similar to above, we choose ﬁve diﬀerent cases: For training, we use monotonic loading and cyclic response for11=0:25 %(50 cycles), 11=0:8 %(3 cycles), and11=1:2 %(3 cycles). As validation, we use 11=0:4 %with 25 cycles. The number of cycles were chosen such that the hysteresis loops are approximately stable during the last cycle. For the mini-batch updates, only ﬁve (out of 18) steps were used for the monotonic loading, ten cycles for 11=0:25 %and one cycle for 11=0:8 %and11=1:2 %to reduce the computa- tional time slightly. Figure 4 shows the evolution of loss during the training process for the neural network and Chaboche models. For the Chaboche model, the training and validation losses de- creaserapidlybeforeremainingconstantduringeachtrain- ing phase. After the mini-batch training phase (ﬁrst 5000 epochs), an additional slight decrease in the losses occur. Thetraininglossfortheneuralnetworkmodelshows(on average) a continuous reduction in the loss, initially rela- tively rapid before slowing down. The validation loss for the neural network model reduces to an initial minimum within the ﬁrst 50 epochs before increasing again during the following 250 epochs. Thereafter, it slowly decreases following the trend of the training loss in the mini-batch updates. In the full-batch updates, it reduces further be- fore rising toward the end. This increase indicates overﬁt- ting, and the parameters from epoch 9467, which had the lowest validation loss, are used. The material in Meyer and Ahlström (2023) exhibits a stress plateau followed by hardening. When analyzing the split between isotropic and kinematic hardening, a rapid decrease in isotropic yield strength combined with a cor- responding increase in back-stress (Fig. 10 in Meyer and Ahlström (2023)). The trained material behavior approx- imates this behavior as well, seen by the hardening mod- uli,Hiso162 GPa andHkin131 GPa (see Table 2), which approximately balance each other, as seen in Figure 5. After the initial plateau phase, where the model shows too high stresses, the ﬁt is almost perfect for the remaining loading. The load case of 11=0:25 %, Figure 6, is hard to ﬁt because the almost fully elastic behavior in the initial cycle. Due to isotropic softening, the plastic strain ampli- tudeincreasesinsubsequentcycles, makingthemodelvery sensitive to the initial behavior. After training, this exper- iment has the highest relative loss of all cases, including the validation case. Despite these challenges, the model can obtain a reasonable loop shape in the 50th cycle. In the training case with 11=1:20 %, Figure 7, the overall 0 2500 5000 7500 1000010−310−210−1 EpochLossTraining Validation Neural net. Neural net. Chaboche ChabocheFigure 4: Loss function evolution for model training on experimental data 0 1 2 3 4 5 6 702004006008001000 /epsilon111[%]σ11[MPa] Experiment Simulation Figure 5: Response of the neural network model for monotonic load- ing (training) −0.2−0.1 0 0.1 0.2−500−2500250500 /epsilon111[%]σ11[MPa]Experiment Simulation Figure 6: Response of the neural network model in the ﬁrst (solid) and 50th cycle (dashed) for =0:25 %(training) 9 loop shapes ﬁt very well with experimental results. Even so, there is a slight waviness compared to the experimental curves. For the validation case, in Figure 8, there are some dis- crepanciesbetweenthemodelresponseandtheexperiment in the second half-cycle. The model shows a more distinct elastic-to-plastic transition compared to the experiments. However, in the 25th cycle, the stress-strain curves are very similar, indicating that the model works quite well for the validation case. 4.2.2. Equation discovery The results in the previous section show that the neural network model can simulate the experimental results rea- sonably well. While this model could be used directly, it contains black-box elements for the evolution laws. Fur- thermore, it is less eﬃcient than a regular material model due to the embedded networks. Instead, we identify sparse symbolic approximations of the network outputs. For 0-3 terms, we use a direct binomial search. This method ﬁnds the best solution since it tests all possible combinations. To be able to ﬁnd expressions with more than three terms, we apply the selection method described in Section 3.2.3. We then run a binomial search on the selected candidate functions to identify up to 6 terms. As a ﬁnal comparison, we also perform a direct sparse regression with a varying penalty factor to obtain expressions with diﬀerent num- bers of terms. The results in Figure 9 show that i) the direct binomial search always ﬁnds the lowest loss for the given numberof terms, and ii)thatlossesforsparse regres- sion followed by the binomial search are only moderately higher. In most cases, the sparse regression followed by a binomial search leads to lower losses than simply increas- ing the penalty parameter until solutions with few enough terms are found. The variability in the sparse regression results highlights the issues with local minima discussed in Flaschel et al. (2021). 4.2.3. Discovered model Based on the results in Figure 9, we evaluate three mod- els closely: i) three terms from direct binomial search, ii) four and iii) six terms from sparse regression followed by a binomial search. Figure 10 shows the training evolu- tion for the two ﬁrst cases compared to the neural net- work model training. Starting with the exact values for the identiﬁed parameters would yield a faster convergence, but all models start with a random initial guess to ensure a fair comparison. Using three identiﬁed terms for each net- work output slightly improves the training loss compared to the Chaboche model and greatly reduces the valida- tion loss (cf. Figure 4). With four terms, the discovered model obtains training and validation losses comparable to the neural network model. Using the identiﬁed model with six terms gives a marginal improvement: The train- ing loss goes from 0:09 %to0:08 %, and the validation loss from 0:08 %to0:07 %. Such marginal improvements do −1.2−0.8−0.4 0 0.4 0.8 1.2−600−400−2000200400600 /epsilon111[%]σ11[MPa]Exp SimFigure 7: Response of the neural network model in the two ﬁrst cycles for=1:2 %(training) −0.4 −0.2 0 0.2 0.4−600−400−2000200400600 /epsilon111[%]σ11[MPa]Experiment Simulation Figure 8: Response of the neural network model in the ﬁrst (solid) and 25th cycle (dashed) for =0:40 %(validation) 0 2 4 6 8 10 1210−310−210−1100101 Number of termsScaled lossNN iso NN k,ν NN k,β Figure 9: Scaled loss for equation discovery using direct binomial search (markers), sparse regression followed by binomial (solid lines), and sparse regression with increasing penalty (dashed). 10 not justify the increased model complexity, and Figure 10 excludes these results for clarity. As four terms seem required for an accurate model, the remaining discussion in this section considers this case, for which the evolution equation for isotropic hardening is _k=_" 1 + Y0+D piso;1+piso;2:p +piso;3h: i+piso;4h: i2E# (37) wherepiso;iare new material parameters. The last term within the outer brackets corresponds to NNisoin Equa- tion (14). Positive network output values, NNiso0, used in the sparse regression, do not guarantee that the approx- imation of NNisois always positive. The expression cor- responding to NNisois, therefore, wrapped in a Macaulay bracket. In addition to ensuring thermodynamic consis- tency, this addition improved the predictive ability of the model of the yield strength evolution in Section 4.2.4. For kinematic hardening, the NNk;andNNk; parts have four terms each, with the complete evolution equa- tion, _b=_" + [: ] pk;1+pk;2jj jj2 +pk;3jj jj+pk;43p h: i  + pk; 1[ :]jj jj2+pk; 4q :p p 2=3 +pk; 3 :p p 2+pk; 2h :i # (38) using the same approach as in Equation (37) to ensure a positive dissipation. Table 2 shows that all models have similar initial yield stress, Y0, and viscoplastic parameters (tandn). While the hardening moduli diﬀer between the models, especially for the Chaboche model, the initial hardening modulus, Hiso+Hkinis within 15 GPafor all models. As expected from the lower training loss, the identiﬁed model ﬁts the monotonic stress-strain curve in Figure 11 Table 2: Trained parameter values for the neural network model (NN), the identiﬁed model with 4 terms (Id.), and the Chaboche model (Ch.). NN Id. Ch. Unit Y0495.0 487.9 486.8 MPa Hiso-162.2 -159.9 -70.02 GPa Hkin131.4 114.3 33.51 GPa t51.87 55.26 71.31 ms n6.885 7.473 6.591 - 0 5000 10000 1500010−310−210−1 EpochLossTraining Validation Neural net. Neural net. 3 terms 3 terms 4 terms 4 termsFigure 10: Loss function evolution for training of identiﬁed models 0 1 2 3 4 5 6 702004006008001000 /epsilon111[%]σ11[MPa] Experiment Chaboche Identified Figure 11: Response of the Chaboche and identiﬁed model for mono- tonic loading (training) −0.2−0.1 0 0.1 0.2−500−2500250500 /epsilon111[%]σ11[MPa]Experiment Chaboche Identified Figure 12: Response of the Chaboche and identiﬁed model in the ﬁrst (solid) and 50th cycle (dashed) for =0:25 %(training) 11 much better than the Chaboche model. However, com- pared to the neural network model in Figure 5, the stress increase deviates more from the experiment in the last two cycles. In the case of 11=0:25 %(Figures 6 and 12), all three models behave practically fully elastic in the ﬁrst cycle. In the 50th cycle, however, the neural network and identiﬁed models give a smoother elastic-to-plastic transi- tion, much closer to the observed experimental behavior than the Chaboche model. For the case of 11=1:20 %in Figure 13, the identi- ﬁed is more accurate than the Chaboche model during the elastic-to-plastic transition and towards the end of each cycle. For this case, the neural network model, Figure 7, exhibited a slight waviness, as previously discussed. This waviness is much more evident for the identiﬁed model, the causes for which are discussed in Section 5.3. In the validation case with 11=0:40 %(Figure 14), the identiﬁed model performs similar to the neural net- work model, with a more distinct elastic-to-plastic transi- tion than the experiment in the second half-cycle, but with an excellent agreement in the 25th cycle. The Chaboche model cannot capture this smoother loop shape and pre- dictsanalmostbi-linearstress-straincurveforbothcycles. 4.2.4. Yield evolution Meyer and Ahlström (2023) measured the yield strength explicitly by probing. Due to the smooth elastic-to-plastic transition for the material, it is generally diﬃcult to accu- rately detect the yield point directly from the stress-strain curve. Consequently, this is also diﬃcult for the mate- rial model to learn based solely on stress-strain data. The primary motivation for the study in Meyer and Ahlström (2023) was to investigate the hypothesis that the isotropic hardening only depends on the accumulated plasticity, . Figure 15 shows that the Chaboche model (circular mark- ers), with Voce hardening, follows the expected gray line, y=Y0+1[1exp(Hiso=1)], for all diﬀerent load cases. However, the experimental results show that the yield strength evolution depends on the loading case. The neural network and the identiﬁed models have simi- lar evolution behavior in all load cases. For 11=0:25 %, the material softens initially before obtaining a constant yield stress. This evolution is qualitatively the same as the experimental results. All models underpredict the plastic- ity in the ﬁrst cycles for this loading case, causing the curves to shift to the left compared to the experiments. For11=0:4 %, the neural network and identiﬁed models correctly predict a somewhat lower yield strength than for 0:25 %but do not fully capture the continued softening in the ﬁrst cycles down to the level of the higher strain amplitudes. For these amplitudes, however, lower yield strengths are correctly identiﬁed, as well as the con- stant level. Forthemonotonicloading, theneuralnetworkandiden- tiﬁedmodelscapturetheinitialsofteningfollowedbyhard- ening. While the models exaggerate this eﬀect compared to the experiments, they could still qualitatively predict −1.2−0.8−0.4 0 0.4 0.8 1.2−600−400−2000200400600 /epsilon111[%]σ11[MPa]Experiment Chaboche IdentifiedFigure 13: Response of the Chaboche and identiﬁed model for the ﬁrst two cycles with =1:2 %(training) −0.4 −0.2 0 0.2 0.4−600−400−2000200400600 /epsilon111[%]σ11[MPa]Experiment Chaboche Identified Figure 14: Response of the Chaboche and identiﬁed model in the ﬁrst (solid) and 25th cycle (dashed) for =0:40 %(validation) 10−410−310−210−11001010100200300400500600 λσy[MPa]/epsilon1=±0.25 % /epsilon1=±0.40 % /epsilon1=±0.80 % /epsilon1=±1.20 % monotonic Figure 15: Yield evolution at cycle peaks for diﬀerent loading types versusaccumulatedplasticity, . Thediﬀerentlinestylesindicateex- periments (solid lines), Chaboche model (markers), neural network (dotted), and identiﬁed model (dashed). The gray background line is the analytical Voce hardening with the Chaboche model’s param- eters. 12 the yield strength behaviors from training on the stress- strain curves. 5. Discussion 5.1. Identiﬁcation of known evolution laws As previously discussed, Figures 2 and 3 show that the neural network model can predict the response for harden- ing laws that the neural network cannot represent exactly. The results in Table 1 showed that the exact evolution laws could be re-identiﬁed by analyzing the trained neural network embedded in the material model. This test shows the potential of the proposed method. 5.2. Capabilities of the neural network model A suﬃciently large feed-forward neural network is a uni- versal function approximator that can describe any func- tion (Hornik et al. (1989)). From this point of view, one would expect that it should be possible to train the neural network model to ﬁt the experimental data exactly, which the results show that it does not do. There are many possible explanations, some of which are discussed in the following paragraphs. The network layout described in Section 2.3 is not a universal function approximator. That is not only due to an insuﬃcient size but also because the last activation function prevents negative values. This constraint ensures fulﬁllment of the 2nd law of thermodynamics and is essen- tial to ensure the physical correctness of the model. We have observed that it improves the model’s predictive abil- ities, even though it increases the training loss slightly, as expected. Using larger networks did not signiﬁcantly improve the results, possibly due to the increased complexity of the non-convex optimization problem. However, this result may also indicate deﬁciencies in the model structure. For example, the standard model assumption of a well- separated elastic and plastic region is only valid for a ho- mogeneous material. Actual test samples consist of a poly- crystalline material in which diﬀerent grains yield at dif- ferent stress levels, giving a smooth elastic-to-plastic tran- sition upon load reversal. Modeling frameworks for rate- independent models with smooth elastic-to-plastic transi- tions were already proposed in the 90s, by e.g. Lubliner et al. (1993). However, the smooth transition also occurs during reloading, see e.g. Lubliner et al. (1993); Rubin and Forest(2020), whichisnotthecasefortheresultsinFigure 5. This transition remains a challenge, even with the very general evolution laws represented by the neural network, highlighting the need for alternative modeling frameworks. 5.3. Interpretation of the identiﬁed behavior Overall, the identiﬁed evolution equation results in a reasonable behavior that is qualitatively similar to the ex- perimental results. The main exception is the waviness observed in Figure 13 for 11=1:20 %. Figure 16 shows /epsilon111 I1=κ I5=ν:/epsilon1p σ11 I3=ν:β 0 9601 1st cycleNormalized values 0 962nd cycle time in cycle [s]Figure 16: Invariant evolution the time evolution of three diﬀerent invariants, the axial strain, and the axial stress for the two ﬁrst cycles. All val- ues are normalized by the maximum absolute value to vi- sualize the evolution of quantities of diﬀerent magnitudes. For comparison, the yield strengths shown in Figure 15 are given byY0+at the end of each cycle. In Figure 16, I1=decreases initially. After the load reverses, it increases rapidly before decreasing again. The remaining cycles follow this cyclic pattern of ini- tial hardening followed by softening. The neural network model gives the same behavior. For the identiﬁed law, the termh: iis responsible for the sudden harden- ing upon load reversal, which makes it possible to model the smooth elastic-to-plastic transition. After the peak of isotropic hardening, a waviness initiates. Many identiﬁed terms contain Macaulay brackets, resulting in discontinu- ousderivativesatzero, contributingtothecontinuedwavi- ness. However, in the remaining experiments and even in other parts of the 11=1:20 %experiment, the behavior is smooth. While the Macaulay terms introduce some challenges in the model, they diﬀerentiate between loading cases. In particular, the diﬀerence in isotropic hardening for plas- tic loading and unloading in Figure 16 seems to be the main factor in predicting the correct qualitative evolu- tion of the yield strength in Figure 15. Very few physi- cal processes show abrupt changes in behavior: While it is reasonable to expect a diﬀerent behavior directly upon load reversal, abrupt changes are unexpected when :p changes sign during monotonic loading. Smoother changes occur by considering, e.g., quadratic terms as :p 2. On the other hand, such terms will grow fast as pin- creases, leading to worse extrapolation. Nonlinear terms with sigmoidal shapes are attractive here (cf. e.g. the tanh activation functions used in the neural network). While 13 those functions are unsuitable for sparse regression, they are compatible with symbolic regression. Udrescu and Tegmark (2020) outperformed all available symbolic re- gression software for the physics problems in “Feynman Lectures on Physics” by ﬁrst training a neural network on the data and then applying symbolic regression on the trained network. That approach aligns with the work pre- sented herein and provides interesting future extensions. 5.4. Implications and opportunities for future studies Whilethepreviousparagraphmentionedpossibilitiesfor improved equation discovery, the present section concerns how to improve the neural network model. The presented model starts from a basic plasticity model with isotropic and kinematic hardening. It is well known that plastic loading can result in both an anisotropic yield surface and stiﬀnessdegradation. Theseeﬀectsmotivateextendingthe presentmodelbyintroducinganevolvingyieldsurfaceand damage evolution. While such extensions are straightfor- ward, suﬃcient suitable experimental training data is re- quired. The suitability of training data is also an important con- sideration when designing the experiments. The predicted yield strength evolution in Figure 16 was unexpected for the authors. This ﬁnding oﬀers essential information for designing experiments to understand yield evolution bet- ter. For example, varying the strain amplitude would pro- vide more valuable data as the cyclic responses stabilize after only a few cycles. 6. Contributions In this paper, we present an approach that embeds a neural network within a conventional material modeling framework. Key advantages of our approach include the intrinsic satisfaction of thermodynamics, frame invariance, and the ability to employ implicit time integration. More- over, we could discover analytical expressions for the evo- lution laws by learning from the trained neural network using both artiﬁcial and experimental data. The iden- tiﬁed laws reveal an interaction between kinematic and isotropic hardening. These new interaction terms enabled the prediction of experimentally observed behaviors of yield strength evolution for various loading cases. 7. Acknowledgements The work in this paper has made use of several open sourcesoftwarepackageswrittenin theJuliaprogramming language (Bezanson et al. (2017)). The automatic dif- ferentiation used ForwardDiff.jl (Revels et al. (2016)), and the implementation relied heavily on the Tensors.jl package (Carlsson and Ekre (2019)). The training used optimizers implemented in Flux.jl , Innes (2018).A. Sensitivity for parameter identiﬁcation Because neural networks have many parameters, a gradient-based optimization algorithm is desired. It re- quires the sensitivity of the loss, Lk(p), in Equation (32), with respect to the material parameters p. AsLkis an explicit function of the simulated time history of stresses, it suﬃces to calculate the derivative of stress wrt. material parameters. Westartbybrieﬂyreviewingtheimplementa- tion of an elasto-plastic material model to make it easier to explain the various functions required to obtain the sought sensitivity. A.1. Material response A general strain-driven, history-dependent, material model, can be considered as two functions1, =m(p;;ns;t) (39a) s=ms(p;;ns;t) (39b) wheresdenotes the hidden material state and tis the time increment. The notationnimplies the value of in the previous time step. For the models considered in this study, the state is given by s= [p;; ]. A standard plasticity model implementation (excluding tangent stiﬀness calculation for brevity) consists of the fol- lowing steps.   1function material_response(p, ,ns,t) 2_trial = elastic_response(p, ,ns) 3 _trial = yield_criterion(p, _trial,ns) 4 if_trial <= 0 # Elastic response 5 return_trial,ns 6 else # Plastic response 7 # Find X such that R(X, p, ,ns,t) = 0 8 X = solve_nonlinear(R, p, ,ns,t) 9 , s = plastic_response(X, p, ,ns,t) 10 return, s 11 end 12end  The trial stress, trial=E: np , on Line 2, is the stress assuming that the current load stress is elastic. If the trial state is inside the current yield surface, 0, the response is elastic, and the trial stress is the correct solution (Line 5). Otherwise, the response is plastic, and thestatevariablesandstressmustbecalculatedbysolving a nonlinear equation system, R(X;p;;ns;t)(cf. Box 7.1 in Neto et al. (2008) for a simpliﬁed example). Given the solution, X(p;;ns;t), which is an implicit function via the nonlinear iterations, the updated state variables and stress can be calculated (Line 9). To later calculate the sensitivity of the loss, dLk=dp, we will need the 6 sensitivities, dmy dz;y2f;sg z2f;ns;pg(40) 1Although in practice calculated together in a single function. 14 The derivatives are trivial and eﬃcient to compute for elastic steps with automatic diﬀerentiation. In the case of a plastic response, however, the implicit function X(p;;ns;t)poses a problem. Therefore, we extend a standard approach for calculating the material tangent stiﬀness. Speciﬁcally, the residual function will always be zero for any input p,, orns, i.e. dR dz= 0 =@R @XdX dz+@R @z(41) With the function plastic_response on Line 9 denoted as^my(X;p;;ns;t), the sought derivatives are dmy dz=d ^my dXdX dz+@^my @z =d ^my dX@R @X1@R @z+@^my @z(42) for plastic loading ( trial>0). The functions Rand ^myare explicit (do not contain any iterations), and their derivativescanbeeﬃcientlycalculatedwithautomaticdif- ferentiation. A.2. Stress state iterations In most experiments, not all stress and strain compo- nents are measured. For the uniaxial and biaxial experi- ments considered in this study, one normal and one shear stress component and the corresponding strain compo- nents are available in the results. The two other normal stresses are zero, while the remaining shear strains are as- sumed to be zero. The corresponding normal strains are unknown, and so are the corresponding shear stresses (al- though these are typically zero). The measured strains are the simulation inputs, and the measured stresses are the outputs. In addition, we have the constraint that the other assumed values (normal stresses and shear strains) must be zero. Hence, from the modeling point of view, we have the known stress components, ~, and the correspond- ing unknown strain components, ~. In this paper, these components are22and33. Additionally, we have the unknown stress components, ^, and corresponding known strain components, ^(components11,12,13, and23 in this study). As the input to the material function is the full strain tensor, an iterative procedure is required to ﬁnd the ~resulting in the prescribed ~. Three diﬀerent outputs are the result of such an iterative procedure, ^=g^(^;ns;p;~) (43a) s=gs(^;ns;p;~) (43b) ~=g~(^;ns;p;~) (43c) where typically ^are the components considered in the loss,Lkin Equation (32). A similar trick as for the material response gives the derivative of the unknown stresses, ^, wrt. the materialparameters: We use that the iterative function will always give the same known stresses, ~, i.e., d~ dp=0=dm~ dp+dm~ d~dg~ dp+dm~ d [ns]d [ngs] dp(44) allowing the derivative of the unknown strains to be cal- culated as dg~ dp=dm~ d~1dm~ dp+dm~ d [ns]d [ngs] dp (45) Given that derivative, the sensitivity of the unknown stresses, ^, can be calculated as dg^ dp=dm^ dp+dm^ d~dg~ dp+dm^ d [ns]d [ngs] dp(46) However, from the last term in Equations (45) and (46), thevalueof d[ngs]=dp, whichstemsfromtheprevioustime step, has yet to be discussed: In the ﬁrst time step, this derivative is zero because the initial state is independent of the material parameters. Assuming that the value in the previous time step is known, the value in the current time step can be calculated as dgs dp=dms dp+dms d~dg~ dp+dms d [ns]d [ngs] dp(47) This result can then be used in the following time step to calculate dg^=dp, and ﬁnally dLk=dp. For performance reasons, it is important to note that for each local iteration in the material model, i.e. solving R(X;p;;ns;t) = 0, only the derivative @R=@Xneeds to be calculated. This derivative is also required without sensitivity calculations. Furthermore, the algorithmic tan- gent stiﬀness, dm=d, must be calculated once per stress stateiteration. Thiscalculationisalsoalwaysneeded. The additional computations for calculating the material sen- sitivities occur only at the end of each time step after the local material model iterations and the stress state itera- tions converge. This possibility reduces the computational costs of obtaining the sensitivity. Declaration of Generative AI and AI-assisted tech- nologies in the writing process During the preparation of this work the authors used ChatGPT in order to improve the spelling, grammar and clarity of the abstract and contribution sections. After us- ing this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.