ON FORWARD AND INVERSE UNCERTAINTY QUANTIFICATION FOR MODELS INVOLVING HYSTERESIS OPERATORS

. Parameters within hysteresis operators modeling real world objects have to be identiﬁed from measurements and are therefore subject to corresponding errors. To investigate the inﬂuence of these errors, the methods of Uncertainty Quantiﬁcation (UQ) are applied. Results of forward UQ for a play operator with a stochastic yield limit are presented. Moreover, inverse UQ is performed to identify the parameters in the weight function in a Prandtl-Ishlinski˘ı operator and the uncertainties of these parameters. intervals, and other Quantities of Interest (QoI) . Moreover, we will also present a brief example of Inverse UQ , i.e. , of using (further) data and measurements, to determine (reduce/adapt) the uncertainty of the parameter(s), i.e. , to determine a (new) random variable taking into account the (new) information, and use the (new) random variable to represent the parameter(s) afterwards.

Forward UQ is presented, for example, in [19] for the Bouc-Wen hysteresis model and in [10] for the Davidenkov model with simplified loading-reloading rules.
Forward and inverse UQ can be found for example in [11] for a model involving the stop operator and in [16] for a model involving a hysteretic multi scale formulation.
Inverse UQ for a Prandtl-Ishlinskiȋ operator of stop type can be found in [13]. The plan of the paper is the following: In Section 2, the play operator with deterministic data is recalled and a method to efficiently compute the values of the play operators for all possible yield limits is introduced. In Section 3, we will consider the play operator if the yield limit is not know exactly such that we represented it by a random variable and perform forward UQ to determine the evolution of the uncertainty. In Section 4, we will recall the generalized Prandtl-Ishlinskiȋ operator with a parameterized weight function as in [4], introduce the corresponding initial loading curve and present equations that are satisfied by this curve. In Section 5, we will consider results of measurements and will try to identify the values of the parameters and extract information about their uncertainties, i.e., we perform inverse UQ. In the first attempt, we extract several approximations of the initial curve from the measurements, use the equations introduced in Section 4 to compute for each approximation a corresponding pair of parameters and investigate the resulting set of pairs. In the second attempt, Bayes' Theorem is applied.

The play operator with deterministic data
An important example for an hysteresis operator is the play operator, see, e.g., [3,8,9,17]. Using a notation as in [3,9,17], we assume that a final time T > 0 and some yield limit r ≥ 0 are given and consider the play operator P r [·, ·] as operator from R × C[0, T ] to C[0, T ] mapping (z 0 , u) ∈ R × C[0, T ] (with z 0 being the initial state, u being the input function and C[0, T ] denoting the space of all continuous functions from [0, T ] to R) to a function P r [z 0 , u] ∈ C[0, T ]. For the function P r [z 0 , u] : [0, T ] → R, we denote by P r [z 0 , u](t) the value at t for any t ∈ [0, T ]. We recall that for u ∈ C[0, T ] being piecewise monotone, i.e., such that there exists 0 = t 0 < t 1 < · · · < t n = T with u being monotone on [t i−1 , t i ] for all i ∈ {1, . . . , n}, it holds that P r [z 0 , u] is piecewise monotone and satisfies P r [z 0 , u](0) = max (u(0) − r, min (u(0) + r, z 0 )) , (2.1a) P r [z 0 , u](t) = max (P r [z 0 , u](t i−1 ), u(t) − r)) , if u is increasing on [t i−1 , t i ], min (P r [z 0 , u](t i−1 ), u(t) + r)) , if u is decreasing on [t i−1 , t i ], (2.1b) for all t ∈ ]t i−1 , t i ] and all i ∈ {1, . . . , n}. An example for an input function and the output of the play operator is shown in Figure 1. To define the function P r [z 0 , u] for general u ∈ C[0, T ], a convergence argument is used, see, e.g., in [8].
In Figure 2, the outputs of the play for one input function and several yield limits are shown. Therein, we observe that even if P r [0, u] and P r [0, u] with r = r coincide on some time interval, they may not coincide in the future.
In the next section, we need to deal with the evolutions of P r [r, u] for all r ≥ 0 simultaneously and are interested in an efficient method to evaluate for t ∈ [0, T ] the function [0, ∞) r → P r [r, u](t). The method in use is inspired by consideration for the Preisach operator and the representation of the Preisach operator by using the play operator and a memory sequence, allowing the evaluation according to Lemma 2.1.
Recalling now II. (2.17) in Proposition II.2.5 from [9] and considering the value of p w,u,t , we see that the assertion is proved.
3. Forward UQ for the play operator with stochastic yield limit Now, we want to consider a situation, wherein the true value of the yield limit is not known exactly, but we think that its value is near to 2. For this situation, we will do forward UQ to determine the influence of this uncertainty.
Hence, we interpret the yield limit r ≥ 0 as value of the random variable R generated from the normal distribution N (2, 0.5 2 ) with mean 2 and standard deviation 0.5 by ignoring (−∞, 0] and rescaling, leading to the following probability density function ρ R of R, see also Figure 3: The mapping [0, ∞) r → P r [w, u](t) is continuous, see e.g., Corollary II.2.2 from [9]. Hence, it follows that the composition of this mapping with R generates a random variable, denoted by For an initial state w ∈ R, an input function u and a time t ∈ [0, T ], we denote by µ w,u,t the probability measure on R for P R [w, u](t) and observe that for the expected value E (P R [w, u](t)) it holds that Hence, we see that the expected value is equal to the output of a Prandtl-Ishlinskiȋ operator with weight function ρ R and an initial state function that is equal to w on [0, ∞).
To compute P R [w, u](t) or the measure µ w,u,t by using (P r [w, u](t)) r≥0 , we consider the representation as in Lemma 2.1.
Let r w,u,t,0 be defined as in (2.2a). Then it follows from Lemma 2.1.a that P r [w, u](t) = w for for all r ≥ r w,u,t,0 . If r w,u,t,0 > 0, we consider (r w,u,t,k ) and p w,u,t as defined in (2.2c)-(2.5). Using is a one-to-one mapping with the derivative (−1) k p w,u,t ∈ {1, −1} such that it holds for r ∈ [r w,u,t,k+1 , r w,u,t,k ] and z ∈ I k that P r [w, u](t) = z if and only if r = z−z w,u,t,k (−1) k pw,u,t . Using this properties, we can show that the probability measure µ w,u,t for the random variable P R [0, u](t) is the sum of a Dirac measure at w weighted by rw,u,t,0 ρ R (r)dr = max{|u(τ )−w|:τ ∈[0,t]} ρ R (r)dr and of a measure with a density on R being equal to 0 if r w,u,t,0 = 0 and being equal to otherwise.

Generalized Prandtl-Ishlinskiȋ operators
Following [6,7,18], we consider generalized Prandtl-Ishlinskiȋ operators that are defined by combining a Prandtl-Ishlinskiȋ operators and a given function g : R → R by mapping an input function u : [0, T ] → [0, T ] to the output of the Prandtl-Ishlinskiȋ considered with the composition g • u used as input function. These kind of operators can be inverted analytically and can model effects like saturation, asymptotic behaviors, etc., see, e.g., discussions in [2,5,18].
In Sections 3 and 5 from [4] the following generalized Prandtl-Ishlinskiȋ operator has been used to model the magnetization of Galfenol for an applied magnetic field H. This operator is defined by mapping Related to a Prandtl-Ishlinskiȋ operator is the initial loading curve, whose value at x ∈ R can be determined by letting the input monotonically increase or decrease linearly from 0 to x on [0, T ] with an initial configuration λ being constant 0, and evaluating the output of the Prandtl-Ishlinskiȋ operator at T .
In applications to plasticity modeled by a Prandtl-Ishlinskiȋ operator this curve can be measured by considering a material with no previous memory, e.g., if the memory was erased for instance by heating above the Curie temperature and cooling again or applying an input function with decreasing oscillations.
The initial loading curve ψ c1,c2 for PI c1,c2 satisfies for all s ∈ R.
Following, [5,6,12], we consider u ∈ C[0, T ] and 0 ≤ t a < t b < t c ≤ T with u being monotone on [t a , t b ] and on [t b , t c ] such that u(t a ) = u(t c ). Recalling (2.1b), we get for all w 0 ∈ R that Considering an input cycle and differences for corresponding measurements for a process that one would like to represent by a Prandtl-Ishlinskiȋ operator, one can derive an approximation for the initial loading curve on [0, |(u(t 2 ) − u(t b ))/2|] from the measurements during [t 1 , t 2 ] by replacing on the right-hand side of (4.5) the values of the operator by corresponding values from the measurements.

Considered situation
In [1], a magnetostrictive Terfenol-actuator is investigated and the hysteresis between the reference signal determining the current generating the magnetic field and the resulting displacement is considered, and the data of the First-Order-Reversal-Curves (FORCs) are used to determine a Preisach operator and a generalized Prandtl-Ishlinskiȋ operator.
We will now use the data measured during the preparation of the FORC-diagram in [1], considering now the measured current (corrente) generated from the reference signal as input data, such that we get Figure 8 from the measurements, quite similar to the FORC diagram in Figure 4a from [1], except for the scale of the input. Now, we will approximate these data by the initial loading curve corresponding to G c1,c2,c3 [H](t), i.e., we consider the differences derived from the function tanh(c 3 I(t)) as input for ψ c1,c1 . This will be used to determine an appropriate value for c 3 and to derive information about c 1 and c 2 and their uncertainty, i.e., for these parameters we will perform inverse UQ.

Identification of c 3
To determine c 3 , we consider the generated approximations of the initial loading curve. For i ∈ {30, · · · , 58}, we have sets of times t i,0 < t i,1 < · · · < t i,Ki such that the current I as function of time has a discrete local maximum at t i,0 and it holds I (t i,0 ) > I (t i,1 ) > · · · > I (t i,Ki ). For i ∈ {1, · · · , 29} it holds that the current I has a discrete local minima at t i,0 and it holds I (t i,0 ) < I (t i,1 ) < · · · < I (t i,Ki ).
Let L denote the relative length change measured by the length sensor as function of time. Then it follows that the data sets (I (t i,k ) , L (t i,k )) Ki k=0 for i = 30, . . . , 58 represent the FORCs. To derive the initial loading curve ζ i,c3 for i ∈ {1, 2, . . . , 58} and for c 3 > 0 by considering the data set (I (t i,k ) , L (t i,k )) Ki k=0 and the transformation x → tanh(c 3 x), we consider (4.5) with ψ c1,c1 replaced by ζ i,c3 , PI c1,c2 [λ 0 , u] replaced by L, u(t) := tanh(c 3 I(t)), t b = t i,0 , and t c = t i,k for all k = 0, . . . , K i , leading to By interpolation, we get a function ζ i,c3 defined on 0, 1 2 |tanh(c 3 I(t i,k )) − tanh(c 3 I(t i,0 ))| . Now, c 3 is determined by requesting that the sum over the squared L 2 -difference between these generated approximations is minimized. We end up with the optimal value c 3 = c tanh = 0.682138. The corresponding (ζ i,c tanh ) 58 i=1 set of initial loading curves is shown in Figure 9.

Identification of c 1,i and c 2,i
For each approximation ζ i,c tanh of the initial loading curve derived from the measurements by using interpolation, we determine (c 1,i , c 2,i ) ∈ (0, 250) × (0.0000001, 25) by minimizing the L 2 -difference between  this approximation for the initial loading curve and the ψ c1,i,c2,i initial loading curve corresponding to the Prandtl-Ishlinskiȋ-operator PI c1,i,c2,i , see an example in Figure 10.
If only the computed parameter pairs for the curves generated for decreasing inputs i.e., the pairs belonging to the FORCs, are considered, and also the pairs belonging to the three shortest definition intervals are ignored, we get the appropriate subset ((c 1,i , c 2,i ))  One could now use this data to start further investigations. As simple example, we consider PI meani(c1,i),meani(c2,i),c3 [u], PI meani(c1,i),meani(c2,i)±stdi(c2,i),c3 [u], and PI meani(c1,i),meani(c2,i)±2 stdi(c2,i),c3 [u], for the initial state λ 0 ≡ 0 and an input function u increasing linearly from 0 to 1.8, decreasing afterwards linearly to 0 and increasing linearly to 1.8 leading to the outputs shown in Figure 11.
As an example for forward UQ, we now assume that the value of c 2 can be represented by a random variable C 2 with with density ρ C2 that is derived from the distribution N mean i (c 2,i ), (std i (c 2,i )) We could consider N m 1 , δ 2 1 and N m 2 , δ 2 2 for some positive constants m 1 , m 2 , δ 1 , δ 2 and could perform some truncation, to derive random variables C 1 and C 2 to be used in UQ to represent an appropriate approximation for c 1 and c 2 . But, this could only be justified if the pairs ((c 1,i , c 2,i )) 58 i=1 or the pairs ((c 1,i , c 2,i )) 55 i=30 could be considered as typical independent samples for N m 1 , δ 2 1 and N m 2 , δ 2 2 . A typical example for the cloud generated by this kind of samples is shown in Figure 13 for N 1, 1 Considering the computed parameter pairs ((c 1,i , c 2,i )) 58 i=1 , see Figure 14, it is obvious that we can not get these set of pairs as samples if we assume that the parameters can be represented by two independent normal distributed random variables, even if the correlation between (c 1,i ) i and (c 2,i ) i is only −0.0737148. Also the 26 data pairs ((c 1,i , c 2,i )) 55 i=30 in the appropriate subset introduced in Section 5.3, see Figure 15, with a corresponding correlation of −0.532375, can not be derived as samples in this way.

Applications of Bayes' Theorem
Now, we will use Bayes' Theorem to derive a probability density on (0, ∞) × (0, ∞) by considering the measured data, aiming to generate C * 1 , C * 2 as discussed above. Let v obs be a vector containing the sets of the non-trivial discrete values for the initial loading curve used to derived their approximations in Sections 5.2 and 5.3, i.e., We consider on (0, ∞) × (0, ∞) the uniform probability density π 0 for (0, 250) × (0.0000001, 25) as priori density. Updating this priori density by following Bayes' Theorem of Inverse Problems, see, e.g., Section 8.1 from [14] and Sections 2.8 and 6.2 from [15], we get the posterior probability density π new on (0, ∞) × (0, ∞) defined by with (c 1 , c 2 )|v obs being the likelihood that the parameter pair has the value (c 1 , c 2 ) if the output v obs is observed. To derive the forward function Ψ c1,c2 : (0, ∞) × (0, ∞) allowing to formulate the likelihood function, we use a method similar to the derivation of (5.1). We consider (4.5) with u(t) := tanh(c 3 I(t)), t b = t i,0 , t b = t i,0 , and t c = t i,k for all k = 0, . . . , K i . Using that L(t i,k ) and L(t i,0 ) are approximations of PI c1,c2 [u](t i,k ) and PI c1,c2 [u](t i,0 ), we end up with the assumption that there is a sample γ = (γ i,k ) we have and that η is a sample for H := and that v obs is a sample for Ψ (c 1 , c 2 ) + H. Hence, we have to consider the likelihood function associated to function mapping (c 1 , c 2 ) ∈ (0, ∞) × (0, ∞) to the random variable Ψ (c 1 , c 2 ) + H.
Hence, we see that this likelihood function will generate a posterior probability density which is numerical equal to 0 except for a small region, such that the 58 computed pairs (c 1,i , c 2,i ) are obviously no typical samples for this probability density.
Recalling the derivation of the likelihood function and [14,15], one realizes that this function is derived from a probability density for a noise disturbed output with using one parameter pair (c 1 , c 2 ) for all possible initial loading curves. Therefore, the Bayes' theorem generated a density on (0, ∞) × (0, ∞) that represents somehow the informations about the position of this one parameter pair that can be extracted by investigating the observed values in v obs and does not contain information for dealing with different parameter pairs for different subsets of the observation, i.e., this application of Bayes' theorem could not generated a density leading to the requested random variable (C * 1 , C * 2 ).

UQ-issues for further research
We would like to find a probability density on (0, ∞) × (0, ∞) such that the 58 computed pairs of parameter values or the appropriate subset of 26 parameter pairs introduced on Section 5.3 are typical samples for for this probability density. A further investigation of Bayes' theorem as discussed for example in [14,15] yields that applying this theorem would require to use these probability densities as parameter instead of the number Figure 17. Normalized likelihood for the standard deviation 0.05 and 0.1 for the noise in the measured length. Left-hand side: 3D plot, right-hand side contour plot with isolines, top: input data range as in Figure 16, bottom: input data range adapted to considered likelihood function.
pairs (c 1 , c 2 ) considered in Section 5.4, and that Bayes' theorem would return a probability density over a set of probability density on (0, ∞) × (0, ∞) as result.
If one should apply now Bayes' theorem or try to determine the probability density on (0, ∞) × (0, ∞) maximizing the likelihood will be subject of further research.

Modeling issues
The considered data are derived from length change over magnetic field data. Hence, to be compatible with the thermodynamical consistent model derived in [4], one would need to replace the generalized Prandtl-Ishlinskiȋ operator G c1,c2,c3 [H](t) by a function evaluation involving also the counter clockwise potential operator of an hysteresis operator appropriate to model the magnetization, e.g., the counter clockwise potential operator to G c1,c2,c3 . Corresponding investigations will be subject of further research.

Conclusion
The output of hysteresis operators depend on parameters, but their values may be not exactly known when modeling real world processes.
These parameter can be considered as random variables and the methods of uncertainty quantification (UQ) allow to investigate the influence of the uncertainty on the output of the model and to determine the random variable representing the information on the parameter.