Sampling Points-Independent Identification of the Fractional Maxwell Model of Viscoelastic Materials Based on Stress Relaxation Experiment Data

Considerable development has been observed in the area of applying fractional-order rheological models to describe the viscoelastic properties of miscellaneous materials in the last few decades together with the increasingly stronger adoption of fractional calculus. The fractional Maxwell model is the best-known non-integer-order rheological model. A weighted least-square approximation problem of the relaxation modulus by the fractional Maxwell model is considered when only the time measurements of the relaxation modulus corrupted by additive noises are accessible for identification. This study was dedicated to the determination of the model, optimal in the sense of the integral square weighted model quality index, which does not depend on the particular sampling points applied in the stress relaxation experiment. It is proved that even when the real description of the material relaxation modulus is entirely unknown, the optimal fractional Maxwell model parameters can be recovered from the relaxation modulus measurements recorded for sampling time points selected randomly according to respective randomization. The identified model is a strongly consistent estimate of the desired optimal model. The exponential convergence rate is demonstrated both by the stochastic convergence analysis and by the numerical studies. A simple scheme for the optimal model identification is given. Numerical studies are presented for the materials described by the short relaxation times of the unimodal Gauss-like relaxation spectrum and the long relaxation times of the Baumgaertel, Schausberger and Winter spectrum. These studies have shown that the appropriate randomization introduced in the selection of sampling points guarantees that the sequence of the optimal fractional Maxwell model parameters asymptotically converge to parameters independent of these sampling points. The robustness of the identified model to the measurement disturbances was demonstrated by analytical analysis and numerical studies.


Introduction
For several decades, fractional-order rheological models have been used to describe, analyze and improve the viscoelastic properties of different materials.In addition to theoretical research dedicated to fractional-order rheological models [1][2][3][4][5], hundreds of studies have been conducted on the applicability of such models for specific materials to describe their mechanical properties.The applicability of such models to the description of different polymers is well known, for example, poly-isobutylene [4], polyurea and PET [6], shape memory polymers [7], amorphous polymers [8] and flax fiber-reinforced polymer [9].Fractional viscoelastic models are also used for modeling laminated glass beams in the pre-crack state under explosive loads [10]; stress relaxation behavior of glassy polymers [11]; description of fiber-reinforced rubber concrete [12]; viscoelastic modeling of modified asphalt mastics [13]; and modeling rate-dependent nonlinear behaviors of rubber polymers [14].The modeling and simulation of viscoelastic foods, for example, food gums [15], carrot root [16], fish burger baking [17], is another field of application of rheological fractional models.Due to the non-integer order of the operations of integration and differentiation, the fractional-order models have improved flexibility and better adjustment to material characteristics, both in the time and frequency domains, compared to those of the classic integer-order models.
Although over the last several decades different fractional differential models have been proposed for modeling the viscoelastic processes in materials, the fractional Maxwell model (FMM) is the best known [4,5].The relaxation modulus of the FMM, described by the product of Mittag-Leffler and inverse power functions, allows for the modeling of a very wide range of stress relaxation processes in materials.Describing the rheological properties of polymers by the FMM [18,19] is well known.However, the FMM was also applied, for example, for computational modeling and analysis on the damping and vibrational behaviors of viscoelastic composite structures [20], viscoelastic flow in a circular pipe [21], effect of temperature on the dynamic properties of mixed surfactant adsorbed layers at the water/hexane interface [22,23] and constitutive equations of the Mn-Cu damping alloy [24].Fractional viscoelasticity described by the Maxwell model turned out to model both exponential and non-exponential relaxation phenomena in real materials.
Different identification methods for the recovery of the parameters of the non-integralorder models, including the FMM, from both static [16,[25][26][27][28] and dynamic [12,[29][30][31] experiments data have been proposed so far.It is known that different identification methods in association with different experiment plans result in different identification data yield models, which may differ [32].Generally, the identification result, i.e., the chosen model, is influenced by the three entries that are necessary for model identification: the set of models from which the best model is chosen, the rule for the optimal model selection and the measurement data obtained in the experiment [32,33].For the selected class of models, here, the set of the fractional Maxwell models, the identified model depends on the identification rule and the experiment data.The model parameters are usually determined by guaranteeing the "best-possible" fit to the measurements.Therefore, parameters of the optimal model are dependent on the measure applied for evaluating the "best" [32].The mean-square approximation error is the predominant selection of the model quality measure, which results in a standard least-squares identification task.For the selected identification index, the model identified is usually dependent, sometimes even very strongly, on the experiment data.This is the case with FMM identification methods known in the literature [12,16,[25][26][27][28][29][30][31].This paper deals with the problem of the FMM identification using measurement data from the stress relaxation test.Therefore, the sampling instants used in the experiment and discrete-time measurements of the relaxation modulus compose the set of the experiment data.To build the optimal fractional Maxwell model whose parameters do not depend on sampling instants applied in the stress relaxation test is the aim of this paper.
In the previous paper [33], the problem of the least-squares approximation of the relaxation modulus has been considered for an assumed wide class of relaxation modulus models.Models being continuous, differentiable and Lipschitz continuous with respect to the parameters have been assumed.The main results in [33] refer to the models that are determined asymptotically, when the number of measurements tend to infinity.Whenever some applicability conditions concerning the chosen class of models are satisfied, the asymptotically optimal FMM parameters can be determined using the measurement data obtained for sampling instants selected randomly due to the appropriate randomization, even when the true relaxation modulus description is completely unknown.For the exponential Maxwell and the exponential stretched Kohlrausch-Williams-Watts models, the applicability conditions are satisfied [33].It should be noted that the concept of identification being measurement point-independent comes from the Ljung paper [34] and the paper of [35], in which the optimal identification problems for dynamic and static systems have been considered.
In this paper, the concept of introducing an appropriate randomization for the selection of sampling instants at which the measurements of the relaxation modulus are recorded is applied for the fractional Maxwell model identification.Following [33], the real material description is completely unknown and only the measurement data of the relaxation modulus are available for model identification.Identification consists of determining the FMM that solves the problem of an optimal least-squares approximation of a real relaxation modulus.The complicated form of the relaxation modulus of the FMM (the product of Mittag-Leffler and inverse power functions) implies that the applicability of the sampling points-independent identification for FMM identification is not obvious.It is known that the relaxation modulus of the FMM is continuous and differentiable with respect to its four parameters [36].However, the satisfaction of the Lipschitz continuous property with the bounded Lipschitz constant is proved in this paper for the first time, to guarantee the applicability of the experiment randomization concept.
A complete identification scheme leading to the strongly consistent estimate of the optimal model was specified.Assuming that the measurements are corrupted by additive disturbances, the stochastic-type analysis of the model convergence was carried out, and the exponential rate of convergence was demonstrated both analytically and by numerical studies.For materials described by the unimodal Gauss-like spectrum of relaxation used to describe the rheological properties of the materials [37][38][39] and by the Baumgaertel, Schausberger and Winter (BSW) spectrum [40,41] successfully applied for modeling the polymers [42,43], based on the simulation experiments, both the asymptotic properties and noise robustness of the algorithm were numerically studied.To improve the clarity of this article, the proof of the new FMM Lipschitz property is moved to Appendix A. The tables with the results of the numerical studies are given in Appendix B.

Material
A linear viscoelastic material subjected to small deformations for which the uniaxial, non-aging and isotropic stress-strain equation is given by a Boltzmann superposition equation [44] is considered, where ς(t) and ε(t) are, respectively, the stress and strain and G(t) denotes the linear relaxation modulus.By Equation (1), the stress ς(t) at time t depends on the earlier history of the strain rate described by the first-order derivative dε(τ) dτ via the kernel given by the relaxation modulus G(t).
The modulus G(t) is the stress induced in the material described by constitutive Equation (1) by the unit step strain ε(t) imposed.It is assumed for the studied material that the mathematical description of the modulus G(t) is completely unknown.However, the real relaxation modulus G(t) is accessible by measurement with a certain accuracy for an arbitrary time t ∈ T .Here, T = [t 0 , T] with the initial time t 0 > 0 and T ≤ ∞.
We make the following assumption [33]:

Fractional Maxwell Model
Constitutive equation of the fractional order Maxwell model is as follows [2,4,45]: where G e denotes the elastic modulus, τ r means the relaxation time, α and β are non-integer positive orders of fractional derivatives of the strain ε(t) and stress ς(t), respectively.In this paper, d α dt α f (x) = D α t f (x) means the fractional derivative operator in the sense of Caputo's of a function f (x) of non-integer-order α with respect to variable t and with a starting point at t = 0, which is defined by [1,4] where n − 1 < α < n and Γ(n) is Euler's gamma function [1] (Equation (A.1.1)).
The FMM (2) can be considered as a generalization of the classic viscoelastic Maxwell model being the series connection of the ideal spring with a dashpot (see Figure 1a) described by a differential equation of the first order [44,46]: with the elastic modulus G e of the spring, the relaxation time τ r = η/G e , where η means the viscosity of the dashpot.
The FMM (2) can be considered as a generalization of the classic viscoelastic Maxwell model being the series connection of the ideal spring with a dashpot (see Figure 1a) described by a differential equation of the first order [44,46]: with the elastic modulus  of the spring, the relaxation time  =   ⁄ , where  means the viscosity of the dashpot.A series connection (see Figure 1c), analogical to the classic Maxwell model, of two elementary fractional Scott-Blair elements ( ,  , ) and ( ,  , ), both described by the fractional differential equation of the general form [2,4,45] with the parameters ( ,  , ) (see Figure 1b), yields the FMM described by Equation (2), where the parameters ( ,  , ) and ( ,  , ) uniquely determine the parameters  and  of the FMM (2); for details, see [16].The four parameters ( ,  , , ) of the FMM (2), compared with only two parameters ( , ), or equivalently ( ,  ) of the classic Maxwell model (3), are important for the improvement in the FMM accuracy and flexibility.
The uniaxial stress response of the FMM (2) imposed by the unit step strain ε(t), i.e., the time-dependent relaxation modulus G(t), for an arbitrary 0 < β < α ≤ 1 is given by the formula [2,4,5]: where E κ,µ (x) is the generalized two-parameter Mittag-Leffler function defined by series being convergent in the whole z-complex plane [1,2]: Further, for the description of the FMM identification task, relaxation modulus model ( 5) is denoted as to emphasize the dependence on a four-element vector of model parameters where the subscript 'M' means the model.For the special case α = β, the FMM (2) reduces to the Scott-Blair model (compare (4)) and the relaxation modulus is described by Let us consider the following set of the FMM admissible model parameters: where β 0 > 0 is an arbitrarily small positive number and the maximal and minimal values of elastic modulus G e and relaxation time τ r follow from the a priori knowledge concerning the material under investigation and are such that G emin > 0 and τ rmin > t 0 .G is a compact subset of the four-dimensional real space R 4 .The properties of the two-parameter Mittag-Leffler function and the model ( 7) have been studied by many authors [1][2][3][4][5].The function E κ,µ (x) (6) is completely monotonic on the negative real axis for 0 < κ ≤ 1 and µ ≥ κ, i.e., the function E κ,µ (−x) is completely monotonic for x > 0, Ref. [4] (Equation (E.32)).Whence, since t 0 > 0, by virtue of (6), for any t ∈ T , and any g ∈ G, we have Let us introduce the function [4] (Equation (E.53)) which, comparing (7) and (13), enables describing the relaxation modulus G M (t, g) (7) in compact form as follows The function e κ,µ (x; λ) (13) is known to play a crucial role in many problems of fractional calculus [4] (p. 372) because it has many excellent and useful properties; some of them were used in this paper.The function e κ,µ (x; λ) is completely monotonic for x > 0 when 0 < κ ≤ µ ≤ 1 whenever the parameter λ > 0 [4] (p. 373) as the product of two completely monotonic functions, which by (14) implies the complete monotonicity of the relaxation modulus model G M (t, g) for t > 0 whenever 0 ≤ β < α ≤ 1.This means, in particular, that for t > 0 and g ∈ G, such that 0 < β < α ≤ 1, the positive definite model G M (t, g) (7) monotonically decreases with increasing t > 0. Therefore, for any t > 0 and any g ∈ G, such that 0 < β < α ≤ 1, in view of ( 12)-( 14), we have where m 0 is defined below by the sequence of inequalities valid for any t ∈ T and any where t 0 > 0.
For the case α = β, the relaxation modulus G M (t, g) (10) is also a completely monotonic function of the time for t > 0, which in view of ( 16) is uniformly bounded for t ∈ T and g ∈ G by G emax m 0 /2.
Therefore, there exists a positive constant i.e., the modulus G M (t, g) is uniformly bounded on the set T × G.
The Lipschitz continuity of the model G M (t, g) with respect to parameter g, which is not obvious, in particular, with respect to non-integer orders of fractional derivatives, is fundamental to guarantee the convergence of the optimal models for the applied here experiment randomization.Therefore, before the identification concept and the respective algorithm are presented, the Lipschitz property of the mapping G M (t, g) (7) will be proved, as summarized in the following theorems.The quite tedious proofs are moved into Appendix A.1.

Lipschitz Continuity of FMM with Respect to Model Parameters
Due to the relation between the parameter α and β, let us consider two cases separately when (a) β < α and (b) β = α.Therefore, the set of admissible model parameters G (11) is decomposed on two disjoint subsets: and in which the relaxation modulus G M (t, g) is described by the formulas ( 7) and (10), respec- tively.The bounded set G 1 is non-closed, i.e., the compactness property of the set G (11) is lost here, while G 2 is compact.
The following spectral representation derived in [47] which results from the known spectral representation of the two-parameter Mittag-Leffler function [1] (Theorem 4.18, Equations (4.7.17) and (4.7.15)) and is valid for 0 < β < α ≤ 1, will be used for g ∈ G 1 .Applying the differential approach in Appendix A.1, the next result is proved.
Theorem 1.Let G 1 defined by (19) be the set of the fractional Maxwell model admissible parameters.Then, the relaxation modulus G M (t, g) (7) of the FMM (2) is continuous and differentiable with respect to g (8) for any time t ∈ T and sup t∈T ,g∈G 1 where ∇ g G M (t, g) denotes the gradient of the function G M (t, g) with respect to the vector g; here, ∥•∥ 2 is the Euclidean norm in the space R 4 .
The above theorem means, in particular, that for an arbitrary small positive β 0 , the mapping G M : T × G 1 → R defined according to Equation ( 7) is, uniformly with respect to the time t ∈ T , a Lipschitz continuous function of the vector of model parameters g with Lipschitz constant M 2 .
In the case (b) β = α, for the set of model parameters G 2 (20), the FMM (2) is described by the power-law relaxation modulus G M (t, g) (10) and the absolute boundness of the gradient ∇ g G M (t, g) is resolved by the next result proved in Appendix A.2. Theorem 2. Let G 2 , defined by (20), be the set of the fractional Maxwell model admissible parameters with equal parametersα and β.Then, the relaxation modulus G M (t, g) (10) of the model ( 9) is continuous and differentiable with respect to g (8) for any time t ∈ T and sup t∈T ,g∈G 2 From the proofs of the above theorems, especially from the nonnegative definiteness of the derivatives (A6) and the two last elements of the gradient ∇ g G M (t, g) (A53), the following property is derived.Property 1.Let G defined by (11) be the set of the FMM (2) admissible parameters.Then, for any fixed time t ∈ T , the relaxation modulus G M (t, g) described by (7) or (10) monotonically increases with increasing parameters G e and τ r and other parameters being fixed, i.e., the greater parameters G e and τ r are, the greater the relaxation modulus G M (t, g) is for the given t ∈ T .

Relaxation Modulus Measurements
Following [33,35], let T 1 , . . ., T N be independent random variables with a common probability density function ρ(t); T is the support of ρ(t).Then, let G i = G(T i ) be the related relaxation modulus of the material for i = 1, . . ., N. Let G i denote their measurements corrupted by additive noise Z i , i.e., G i = G i + Z i , recorded in the stress relaxation experiment [44,46,48].
The two assumptions concerning the measurement noises are taken (compare Assumptions 5 and 6 in [33]) as follows: Assumption 2. The measurement noise {Z i } is a time-independent, i.e., independent of the variables {T i }, sequence of independent identically distributed (i.i.d.) random variables with zero mean E[Z i ] = 0 and a common finite variance E Z 2 i = σ 2 < ∞.Assumption 3. The measurement noises Z i are bounded by δ, i.e., |Z i | ≤ δ < ∞ for i = 1, . . ., N.
Both the above assumptions and Assumption 1, concerning the real relaxation modulus, are natural in the context of the relaxation modulus identification [33].

Identification Problem
FMM identification involves selecting from a given class of models defined by (7) and (10) the model that best fits the measurement data.Suppose an identification experiment resulted in a set of measurements G(T i ) = G(T i ) + Z i at the sampling times T i ≥ t 0 > 0, i = 1, . . ., N. The mean-squares index is taken as a measure of the FMM model accuracy.Here, the lower index denotes the number of measurements.Then, the problem of the optimal model identification consists of the solution of the minimization task where g * N is the optimal model parameter.Since, due to the continuity of the model G M (t, g) with respect to the parameter g, the index Q N (g) is a continuous function of g and the set of admissible parameters G (11) is compact, the existence of the solution to the optimization problem (25) immediately results from the Weierstrass theorem about the extreme of continuous function on the compact set [49].Since the minimum g * N can be not unique, let G * N denote the set of vectors g * N that solve the optimization task (25).The parameters g * N of the identified relaxation modulus model G M t, g * N are dependent on the measurement data, in particular, on the sampling instants T i .To make the model independent of specific sampling instants T i , the optimal sampling points-independent approximation problem is stated in the following subsection.

The Optimal FMM
Let us consider the following problem of determining such an FMM that minimizes the global approximation error: where the selected weight function, such that 0 ≤ ρ(t) ≤ M 0 < ∞, is a density on the set T , i.e., T ρ(t)dt = 1.The integral (26) is absolutely integrable, uniformly on G, both for the bounded or unbounded domain T as the product of a function [G(t) − G M (t, g)] 2 , in view of (18)  bounded uniformly for (t, g) ∈ T × G, and absolutely integrable function ρ(t).Therefore, the integral ( 26) is well defined for any g ∈ G.
The problem of the optimal approximation of the real modulus G(t) within the class of the fractional Maxwell models relies on determining the parameter g * that minimizes Q(g) over the set G, i.e., in solving optimization task Due to continuity of G M (t, g) with respect to the vector g, the index Q(g) ( 5) is a continuous function of g, and thus, the existence of the solution g * follows from the previously mentioned Weierstrass theorem concerning the extreme of continuous function on the compact set.Let the set of model parameters g * solving (27) be denoted by G * .Any parameter g * ∈ G * does not depend on the specific time instants applied in the experiment.

Results and Discussion
In this section, the analysis of the asymptotic properties of the identified fractional Maxwell model, when the number of measurements tend to infinity, is conducted.The rate of the convergence of this model to the optimal FMM, which does not depend on the experiment data, is studied.The resulting identification algorithm is outlined.Next, the analytically proven properties of the identification method are verified by numerical simulations and studies.Two example materials are simulated.In the first, the "real" material is described by a unimodal Gauss-like relaxation spectrum [37][38][39] with short relaxation times and the Baumgaertel, Schausberger and Winter (BSW) spectrum [40,41] with long relaxation times.Such models are used to describe the rheological properties of various materials, especially polymers and biopolymers.Based on the noise-corrupted data from the simulated randomized stress relaxation experiment, the optimal FMM models are determined.The asymptotic properties and noise robustness have been studied.

Convergence
The empirical index Q N (g) (24) can be obtained by the replacement of the integral in Q(g) (26) with the finite mean sum of squares, which is clear from a practical point of view.For i = 1, . . ., N, by Assumption 2, the expected value is whence, by (24), the expected value is To investigate the stochastic-type asymptotic properties of the empirical identification task given by ( 25), some properties derived in [35] will be used.Note, that Assumptions A1-A3 from [35], concerning the compactness of the set of model admissible parameters, continuity, differentiability and Lipshitzness of the model are satisfied here.Taken above, Assumption 2 is identical with Assumption A5 in [35], while property ( 18) also means that Assumption A4 from [35] is satisfied, i.e., all the assumptions from [35] hold here.
Property 2. When Assumptions 1 and 2 are satisfied, then where w.p.1 means "with probability one".
By ( 28) and (29), the empirical identification index Q N (g) (24) is arbitrarily close to its expected value, uniformly in g over the set G. In consequence, the model parameter g * N solving empirical identification task (25) can be related to the parameter g * that solves the sampling points-independent minimization task (27).From the uniform in g ∈ G convergence of the index Q N (g) in (29), we conclude immediately the main result of this subsection, c.f., Assertion in [35] or Equation (3.5) in [34].Property 3. Assume that Assumptions 1 and 2 hold, T 1 , . . ., T N are independently and randomly selected from T , each according to the distribution with probability density function ρ(t).If the solutions to the minimization problems ( 25) and ( 27) are unique, then for all t ∈ T .If the minimization problems ( 25) and ( 27) do not have unique solutions, then for any convergent subsequence of the sequence g * N , where and for any t ∈ T and some g * ∈ G * ,the convergence in (31) holds.
The existence of a convergent subsequence of g * N so that the asymptotic property (32) holds results directly from the compactness of G (11).Therefore, under Assumptions 1 and 2, the optimal parameter g * N of the FMM is a strongly consistent estimate of some parameter g * ∈ G * .
Since, by Theorems 1 and 2, the model G M (t, g) is Lipschitz on G uniformly in t ∈ T , then the almost-sure convergence of g * N to g * in (30) implies that, c.f., (Ref.[35]: Remark 2): i.e., that G M t, g * N is a strongly consistent estimate of the optimal FMM G M (t, g * ), uni- formly on T .
Concluding, when Assumptions 1 and 2 are satisfied, the arbitrarily fine approximation of the FMM with the optimal parameter g * can be determined (almost everywhere) as the number of measurements N grow enough, even if the real description of the material modulus is fully unknown.

Exponential Rate of Convergence
Analyzing the convergence in ( 30) and ( 32), the question immediately arises of how fast g * N tends to some g * ∈ G * as N grows large.As in [35], the distance between the model parameters g * N and g * will be evaluated by means of the integral identification index Q(g) (26), i.e., in the sense of the difference Q For this purpose, it will be checked how fast, for a given small ε > 0, the probability P Q(g * ) − Q g * N ≥ ε tends to zero, as N increases.From the well-known Hoeffding's inequality [50], the upper bound of this probability can be derived, analogous to inequality (15) in [35] or inequality (22) in [33] (for details, see Appendix A.1 in [33]): for any ε > 0, where with the constants M and M 1 defined in Assumption 1 and Equation ( 17), respectively, the noises' variance σ 2 and upper bound δ are introduced by Assumptions 2 and 3.
The inequality (34) describes the influence of the number of measurements N and the noises' "strength" on the rate of convergence.For ε being fixed, the bounds for P Q(g * ) − Q g * N ≥ ε decrease exponentially to zero as N increases.The convergence rate is the higher, the lower is M (35).In particular, a quick inspection of (35) shows that for stronger measurement noises, the rate of convergence is reduced.Larger δ and σ yield a greater decrease in the rate.This is as expected, since with large disturbances, the measurements are not very adequate.Simultaneously, the larger M + M 1 , i.e., in view of the estimation (18), the greater the discrepancy between the real modulus and the FMM, the worse the convergence.

Identification Algorithm
In view of the convergence properties (30), (31) the computation of the parameter g * N approximation the parameter g * of the optimal FMM requires the next steps: 1.
Select randomly from the set T the sampling times t 1 , . . ., t N , choosing each t i independently, according to the probability distribution of the density ρ(t) defined given by the weight function in the integral Q(g) (26).

2.
Conduct the stress relaxation experiment [44,46,48], measure and store the measurements G i of the relaxation modulus for the selected time instants t i , i = 1, . . ., N.

3.
Solve the identification optimization task (25) and compute the identified model parameter g * N .4.
Put N = N and g * N = g * N .To extend the set of experiment data, select new N ≫ N.

5.
Repeat Steps 1-3 for a new N, that is, randomly choose new sampling times, conduct the rheological experiment once more for a new sample of the material and determine the next g * N .6.
Examine if ∥g * N − g * N ∥ 2 < ε, where ε is a small positive number, to check if g * N is an adequate approximation of g * .If yes, stop the scheme and take g * N as the approximate value of g * .Otherwise, go again to Step 4.

Remark 1. A less restrictive testing regarding whether |Q
| < ε holds can be used as an alternative for the stopping rule from Step 6.Both types of stopping rules are commonly used in numerical optimization techniques.

Numerical Studies
The results of the numerical studies are concerned with the asymptotic properties of the determined optimal FMM and the influence of the measurement noises on this model.Apart from the theoretical analysis above, these simulation studies make it possible to show the respectability and effectiveness of the method developed for FMM identification.
The "real" material and the FMM model were simulated in Matlab R2023b, The Mathworks, Inc., Natick, MA, USA.Functions MLFFIT2 [58] and MLF [59], provided by Podlubny, were used for the FMM simulation and numerical solution of the optimal identification tasks.

Material I
Consider the material whose relaxation spectrum is described by the unimodal Gausslike distribution: where the parameters are as follows [60]: ϑ = 31520 Pa•s, m = 0.0912s −1 and q = 3.25 × 10 −3 s −2 .The related relaxation modulus is [60] where the complementary error function er f c(x) is given by [4] (Equation (C.2)) Following [47], for numerical simulations, the time interval T = [0, 200] seconds is chosen.Hence, the weighting function in The elements of the optimal parameter vector g * solving the measurement-independent optimization task (27) are given in Table 1.
Table 1.The components α * , β * , G * e and τ * r of the FMM parameter g * solving the optimal identification problem ( 27) and the optimal integral quadratic indices Q(g * ) defined by (27) for the "real" relaxation modulus G(t) (36).The N sampling instants t i for the simulated stress relaxation test were selected randomly according to the uniform distribution on T .A normal distribution with zero mean value and variance σ 2 was applied to the random independent generation of the additive measurement noises {z i }.In the noise robustness analysis, the standard deviations σ = 2, 5, 8 [Pa] were used.In the analysis of the model asymptotic properties, for any σ num- bers of measurements, N ∈ N have been applied, where N = {50; 100; 200; 500; 1000; 2000; 5000; 7000; 10,000; 12,000; 15,000} .

Asymptotic Properties
Then, for every pair (N, σ), the optimal parameter g * N was determined through solving the minimization task (25).The elements of the vectors g * N , the mean square indices Q N g * N and integral Q g * N indices, and the relative percentage errors of the approximation of the measurement-independent parameter g * , defined as are given in Tables A1-A3 for the three standard deviations of the noises.The model approximation error was also estimated via the relative mean error defined as (compare (24)) The optimal model parameters g * N as the functions of the number of measurements N are illustrated by Figure 2 for the noises of σ = 2, 5, 8 [Pa].In any subplot, the values of the related parameters of the sampling points-independent model g * are depicted by horizontal purple lines.The asymptotic properties are also illustrated by Figure 3 juxtaposing the empirical index mean-square index Q N g * N , Equation (24), and the integral quadratic sampling instants-independent index Q g * N , Equation ( 26), as the functions of N with the index Q(g * ), marked with horizontal lines.In Figures 2 and 3, a logarithmic scale is applied for the horizontal axes.These plots confirm the asymptotic properties of the proposed identification algorithm.The convergence of g * N to the parameter g * is directly translated into the convergence of Q g * N into Q(g * ), especially for N ≥ 5000.The values of the index Q N g * N for N = 50, small compared to those for N ≥ 100 (see Tables A1-A3), result from the good fit of the FMM, whose four parameters are optimally selected in problem (25), to only 50 measurement points.For more measurement points, such a good fit is, generally, impossible whenever the real characteristic does not depend on the pre-assumed class of models.A comparison of Figures 2 and 3b with Figure 3a shows that the impact of stronger noises on the values of the empirical index Q N g * N is much stronger than the impact of the noises on the values of the FMM parameter g * N and, consequently, also on the integral index Q g * N , which does not directly depend on the measurement noises.Given Equation ( 28), this property is natural and fully justified.The quality of the real modulus () approximation by the FMM is illustrated Figure 4, where the measurements  ̅ of the real modulus () fitted by the optim model  (,  * ) are plotted for the  = 100 and  = 10,000 measurements and t strongest disturbances;  = 8 [Pa].Although, for the  = 100 measurements, the mo els  (,  * ) and  (,  * ) differ slightly (see small subplot), for the  = 10,0 measurements, they are practically identical, which is confirmed by the values of  (37) equaling 0.52% for  = 100 and equaling only 5.58 × 10 −4 % for  = 10,000 (s Table A3).Even for the strongest noises, the relative errors  (37) of the paramet  * and  * discrepancy is smaller than 0.002% for  ≥ 200.This almost excellent fitti  ; the horizontal purple lines correspond to the optimal integral index Q(g * ) defined in Equation (27).
The quality of the real modulus G(t) approximation by the FMM is illustrated in Figure 4, where the measurements G i of the real modulus G(t) fitted by the optimal model G M t, g * N are plotted for the N = 100 and N = 10, 000 measurements and the strongest disturbances; σ = 8 [Pa].Although, for the N = 100 measurements, the models G M t, g * N and G M (t, g * ) differ slightly (see small subplot), for the N = 10, 000 measurements, they are practically identical, which is confirmed by the values of ERR (37) equaling 0.52% for N = 100 and equaling only 5.58 × 10 −4 % for N = 10, 000 (see Table A3).Even for the strongest noises, the relative errors ERR (37) of the parameters g * and g * N discrepancy is smaller than 0.002% for N ≥ 200.This almost excellent fitting of the experiment data by the model G M t, g * N is confirmed by the values of the relative square model approximation index Q Nrel g * N (38), which for N ≥ 200 and the weakest noises does not exceed 0.015%, while for the strongest noises, it does not exceed 0.28%.

Noise Robustness
To examine the effect of the measurement noises, for every pair (, ), the simulated experiment was repeated  = 50 times.In each experiment repetition, the measurement noises  were generated independently and randomly with a normal distribution, with a zero mean value and variance  .
Having in mind the definition of the index  () (38), for the -element sample, the mean relative relaxation modulus approximation error was determined as follows: for any pair (, ), where the vector of the optimal FMM parameters  , * was computed for j-th experiment repetition,  = 1, … , .
For the true relaxation modulus approximation, the mean optimal integral error was determined for every pair (, ).The distance between the vector  , * and the measurement-independent vector  * for the  element sample was estimated by the mean relative error defined as follows (compare  (37)): The indices  (39) and  (40), as the functions of  and , are depicted in the bar in Figure 5, while the index  (41) is shown in Figure 6.

Noise Robustness
To examine the effect of the measurement noises, for every pair (N, σ), the simulated experiment was repeated n = 50 times.In each experiment repetition, the measurement noises {z i } were generated independently and randomly with a normal distribution, with a zero mean value and variance σ 2 .
Having in mind the definition of the index Q Nrel (g) (38), for the n-element sample, the mean relative relaxation modulus approximation error was determined as follows: for any pair (N, σ), where the vector of the optimal FMM parameters g * N,j was computed for j-th experiment repetition, j = 1, . . ., n.
For the true relaxation modulus approximation, the mean optimal integral error was determined for every pair (N, σ).The distance between the vector g * N,j and the measurement-independent vector g * for the n element sample was estimated by the mean relative error defined as follows (compare ERR (37)): The indices ERRQ Nrela (39) and ERRQ (40), as the functions of N and σ, are depicted in the bar in Figure 5, while the index MERR (41) is shown in Figure 6.
From Figure 5b, it is seen that for N > 2000, the number of measurements do not essentially affect the integral index ERRQ, either for weak or strong noises, while both the empirical index ERRQ Nrel and mean relative error MERR decrease exponentially with the increasing number of measurements, which confirms the analytical analysis performed above.The MERR index is of order 0.55% for N = 100, it does not exceed 10 −3 % for N ≥ 1000 and is smaller than 5 × 10 −5 % even for the strongest disturbances.This, practically, means determining the sampling points-independent parameter g * .The algorithm ensures the very good quality of the measurement approximation even for large noises.The values of the relative relaxation modulus approximation error ERRQ Nrel , which due to the "real" modulus model difference is lower bounded by 3.191 × 10 −4 %, already for N ≥ 100 measurements do not exceed 0.35%, and for N ≥ 1000, fall below 0.028%.The course of the mean integral sampling points-independent index ERRQ (40) as the function of N indicates the asymptotic independence of the model from the sampling points, especially for N ≥ 5000.From Figure 5b, it is seen that for  > 2000, the number of measurements do not essentially affect the integral index , either for weak or strong noises, while both the empirical index  and mean relative error  decrease exponentially with the increasing number of measurements, which confirms the analytical analysis performed above.The  index is of order 0.55% for  = 100, it does not exceed 10 % for  ≥ 1000 and is smaller than 5 × 10 % even for the strongest disturbances.This, practically, means determining the sampling points-independent parameter  * .The algorithm ensures the very good quality of the measurement approximation even for large noises.The values of the relative relaxation modulus approximation error  , which due to the "real" modulus model difference is lower bounded by 3.191 × 10 %, already for  ≥ 100 measurements do not exceed 0.35%, and for  ≥ 1000, fall below 0.028%.The course of the mean integral sampling points-independent index  (40) as the function of  indicates the asymptotic independence of the model from the sampling points, especially for  ≥ 5000.From Figure 5b, it is seen that for  > 2000, the number of measurements do not essentially affect the integral index , either for weak or strong noises, while both the empirical index  and mean relative error  decrease exponentially with the increasing number of measurements, which confirms the analytical analysis performed above.The  index is of order 0.55% for  = 100, it does not exceed 10 % for  ≥ 1000 and is smaller than 5 × 10 % even for the strongest disturbances.This, practically, means determining the sampling points-independent parameter  * .The algorithm ensures the very good quality of the measurement approximation even for large noises.The values of the relative relaxation modulus approximation error  , which due to the "real" modulus model difference is lower bounded by 3.191 × 10 %, already for  ≥ 100 measurements do not exceed 0.35%, and for  ≥ 1000, fall below 0.028%.The course of the mean integral sampling points-independent index  (40) as the function of  indicates the asymptotic independence of the model from the sampling points, especially for  ≥ 5000.

Material II
Consider the material described by the empirical spectrum of relaxation times  introduced by Baumgaertel, Schausberger and Winter [40,41], Figure 6.Dependence of the mean relative error MERR (41) between the optimal parameters g * N and g * of the FMM approximating the "real" relaxation modulus G(t) (36) on the number of measurements N and the noises' standard deviation σ.
Table 2.The elements α * , β * , G * e and τ * r of the FMM parameter g * solving the optimal identification task ( 27) and the optimal integral quadratic indices Q(g * ) defined by (27) for the "real" relaxation modulus G(t) (42), (43).As previously described, in the simulations, the sampling points t i were randomly selected according to the uniform distribution on T .The standard deviations σ = 3, 6, 8 [kPa] of the random normally distributed noises {z i } combined with the number of measure- ments N ∈ N were used for the analysis of the model asymptotic properties.

Asymptotic Properties
For every pair (N, σ), the elements of the optimal model parameter g * N , the empirical and integral Q g * N indices and the relative percentage errors ERR (37) are given in Tables A4-A6 in Appendix B. The dependence of the optimal model parameters g * N on the number of measurements N for the noises of σ = 3, 6, 8 [kPa] are illustrated by Figure 7. Figure 8 illustrates the empirical Q N g * N and integral Q g * N indices as the functions of N; the value of Q(g * ) is marked with purple horizontal lines.These plots confirm the asymptotic properties of the proposed identification algorithm.Figure 8a shows the impact of noises on the values of the empirical index Q N g * N .
Materials 2024, 17, x FOR PEER REVIEW 18     24), (b) the integral quadratic sampling instantsindependent index Q g * N , Equation ( 26), as the functions of the number of measurements N for the noises σ = 2, 5, 8 [Pa]; the horizontal purple lines correspond to the optimal integral index Q(g * ) defined in Equation (27).
The approximation of the real modulus G(t) by the FMM is illustrated in Figure 9, where the measurements G i of the real modulus G(t) along with optimal models G M t, g * N and G M (t, g * ) are plotted for the N = 100 and N = 10, 000 measurements and the noises σ = 8 [Pa].However, for N = 100, the model parameter error ERR = 5.35%, while for N = 10, 000, we have ERR = 0.15%; both for N = 100 and N = 10, 000, the models G M t, g *

Noise Robustness
For every pair (, ), the simulated experiment was repeated  = 50 times.The mean relative relaxation modulus approximation error  (39), the mean optimal integral error  (40) and the mean relative error of the parameter  * approximation  (41) were determined.The indices  and  are depicted in Figure 10 as the functions of  and . Figure 11 illustrates the dependence of the index  on  and .

Noise Robustness
For every pair (N, σ), the simulated experiment was repeated n = 50 times.The mean relative relaxation modulus approximation error ERRQ Nrel (39), the mean optimal integral error ERRQ (40) and the mean relative error of the parameter g * approximation MERR (41) were determined.The indices ERRQ Nrela and ERRQ are depicted in Figure 10 as the functions of N and σ.The mean integral error  for  ≥ 12,000 does not depend essentially on the number of measurements, either for small or large noises (see Figure 10b), while both the empirical index  and mean relative error  decrease exponentially with the increasing number of measurements, the  for  ≥ 7000.For  ≥ 1000, the  index does not exceed 1.01%, for  ≥ 7000 it does not exceed 0.22%, while for  ≥ 10,000, it falls below 0.05% even for the strongest disturbances.The globally optimal   The mean integral error ERRQ for N ≥ 12, 000 does not depend essentially on the number of measurements, either for small or large noises (see Figure 10b), while both the empirical index ERRQ Nrel and mean relative error MERR decrease exponentially with the increasing number of measurements, the MERR for N ≥ 7000.For N ≥ 1000, the MERR index does not exceed 1.01%, for N ≥ 7000 it does not exceed 0.22%, while for N ≥ 10, 000, it falls below 0.05% even for the strongest disturbances.The globally optimal parameter g * was determined.As is seen from Figure 9, the algorithm practically ensures an excellent quality of the relaxation modulus approximation even for the strongest noises.The values of the relative relaxation modulus approximation error ERRQ Nrel , already for the N ≥ 100 measurements, do not exceed 3.3 × 10 −4 % and for N ≥ 1000, fall below 8.3 × 10 −6 %.From the course of the mean integral sampling points-independent index ERRQ (40), as the function of N, we can conclude that the model is practically independent on the sampling instants for N ≥ 12, 000, independently on the measurement noises.The above combined with the close to zero values of ERRQ Nrel means the determining of the globally optimal model with the parameter g * .In conclusion, the courses of both the index ERRQ Nrel (38), and the indices MERR (41) and ERRQ (40) as the functions of N, indicate the asymptotic independence of the model from the sampling points for a sufficiently large number of measurements.

Conclusions
The fractional Maxwell model allows for the modeling of a very wide range of stress relaxation processes in materials.The goal of the FMM identification is, generally, not to achieve a true description of the genuine relaxation modulus, but one that is "optimally accurate" in the assumed sense of the square weighted approximation error and does not depend on the particular sampling instants used in the stress relaxation experiment.The stochastic-type analytical analysis and numerical studies demonstrated that, despite the fact that the real description of the relaxation modulus is completely unknown, an arbitrarily exact approximation of the sampling points-independent optimal FMM can be identified based on the relaxation modulus data sampled randomly, according to respective randomization, when the number of the measurements applied in the experiment appropriately grow large.The four parameters of the approximate FMM are strongly consistent estimates of the parameters of the sampling points-independent model minimizing the integral square approximation error.The resulting identification scheme is simple and useful in application.It requires only the a priori, before the experiment is performed, independent random choice of the time instants at which the relaxation modulus is recorded from the assumed set according to a stationary rule.
Although this article is about modeling the relaxation modulus, the proposed identification scheme can also be successfully applied to the of the fractional-order models of creep compliance using the measurements obtained in the retardation test, whenever the respective set of sampling instants is open to manipulation during experimental data collecting.Therefore, the applicability of the identification asymptotically independent of the time instants used in the rheological experiment, to other fractional-order models determination, in particular, Kelvin-Voight, Zener and anti-Zener models, can be the subject of future research.
To estimate the value of (A7) for t = t 0 , the series representation resulting directly from ( 6) and (A7), is used.The first summand of the series is positive, while the next elements are positive or negative, depending on the index n and the relation between parameters α and β.Since and the argument of the gamma function is such that which in view of the monotonicity of the gamma function implies where x min ∼ = 1.4616321 is the real nonnegative argument at which a minimum of the function Γ(x) occurs [64].In view of (A8) and (A9), having in mind the nonnegative definiteness of , we obtain the following estimation , which, by the inequalities τ rmin ≤ τ r ≤ τ rmax , for t = t 0 , implies the next estimation By the assumption t 0 < τ rmin , the above estimation can be rewritten in compact form as Since, for an arbitrary β 0 ≤ β < α ≤ 1, we have i.e., the derivative for t = t 0 is bounded, uniformly on the set G 1 , where positive parameter m 0 is defined in Equation ( 16).Since the continuous function of the time t is bounded for t = t 0 and for t → ∞ for any fixed g ∈ G 1 , derivative ∂G M (t,g) ∂τ r as a function of the time is bounded both for the bounded and not-bounded set T .However, due to the non-compactness of the set G 1 , from which α = β is excluded, the uniform on T × G 1 boundness of ∂G M (t,g) ∂τ r is not obvious.Therefore, it should be examined if the maximum of (with respect to the time) is bounded, as α → β + .To this end, an alternative to (A6) and (A7), the representation of is derived based on the spectral representation given by Equation ( 21).Differentiation of (21) on both sides with respect to τ r yields where the function by (21), is such that Whence, introducing the notations Equation (A11) can be rewritten as or in a more compact form as a linear combination: where the integrals: The denominator of the fractions r(v, g), r 1 (v, g) and r 2 (v, g) is positive for any v ≥ 0, whenever α − β ̸ = 1, i.e., for any admissible parameter g ∈ G 1 , which satisfies the following inequalities For α → β + , the denominator q(v, g) → [(τ r v) α−β + 1] 2 ≥ 1 .By (A13), (17) and nonnegative definiteness of r(v, g) (A12) on the set R + × G 1 , we have for any (t, g) ∈ T × G 1 , i.e., the function r(v, g)v α−1 e −tv as the function of the variable v is absolutely integrable, uniformly on the set T × G 1 , with the constant m (compare (A2)) given by where the parameter 0 < γ 2 = min In view of Property A2, bearing in mind inequality (17), or (15), to prove the absolute uniform boundness of the derivative ∂G M (t,g) ∂τ r (A16) on the set T × G 1 , it is enough to demonstrate that the integrals I 1 (t, g) (A17) and I 2 (t, g) (A18) are convergent and ab- solutely bounded, uniformly on the set T × G 1 .For this purpose, we express the two integrals as definite integrals of the product of some absolutely integrable function and a bounded function.and according to Property A2, the integral I 2 (t, g) (A18) is convergent for any t ∈ T and any g ∈ G 1 and absolutely bounded uniformly on T × G 1 , with the upper bound G emin γ 2 M 1 resulting from (A20) and (A28).Combining the absolute of the three summands of the right-hand side of (A16), uniform on the set T × G 1 , the respective uniform boundness of ∂G M (t,g) ∂τ r is proved.21) on both sides with respect to β yields where r(v, g) is given by (A12), which using the notations r 1 (v, g) (A14) and r 2 (v, g) (A15) and introducing functions can be rewritten in compact form as a linear combination: of four integrals: To prove the absolute boundness of the derivative ∂G M (t,g) ∂β (A32), uniform on the set T × G 1 , it is enough to demonstrate that the four above integrals are convergent and absolutely bounded, uniformly on T × G 1 .
For any (v, g) ∈ R + × G 1 , the continuous function r 3 (v, g) (A30) satisfies the following inequalities (compare (A14) and (A22)) where m 1 is defined in (A22), i.e., r 3 (v, g) is absolutely bounded, uniformly on R + × G 1 , which combined with the absolutely integrability of the function v α−1 e −tv for any t ∈ T and any g ∈ G 1 , according to Property A2, implies the convergence of the integral I 3 (t, g) (A33) for any t ∈ T and any g ∈ G 1 and its absolute boundness by the upper bound m 1 Γ(β 0 )/γ 1 resulting from and (A24), (A25) and (A37), valid uniformly on the set T × G 1 .
To examine the properties of r 5 (v, t, g) (A43), the asymptotic properties as v → 0 + and v → ∞ are studied.This function is expressed as , where the nominator tends to −∞ and the denominator tends to +∞, as the variable v → 0 + .Therefore, by applying the L'Hospital's rule, in view of α > β > 0, we obtain lim Since (τ r v) α−β ln(τ r v) tends to +∞, while e − 1 2 tv tends to zero, as the variable v tends to infinity, using the L'Hospital's rule double times to the right expression in (A43), we have being valid for any integer n, k > 0 and real x > 0, by putting n = 1 and k = 2 function r 5 (v, t, g) (A43) for any (v, t, g) ∈ R + × T × G 1 can be estimated by For τ rmax v < 1, the right inequality in (A47) implies while for τ rmax v ≥ 1, the middle inequality in (A47) yields where positive m 0 is defined in Equation ( 16).Combining (A48) and (A49), we obtain the inequality valid for any (v, t, g) ∈ R + × T × G 1 , which together with (A44) and (A45) means an absolute boundness of continuous functions r 12 (v, g)r 5 (v, t, g) and r 22 (v, g)r 5 (v, t, g), re- spectively, and in view of the absolute integrability of v α−1 e − 1 2 tv and r(v, g)v α−1 e − 1 2 tv , imply the convergence of the integrals I 5 (t, g) (A39) and I 6 (t, g) (A40) for any t ∈ T and any g ∈ G 1 .The absolute boundness of I 5 (t, g) and I 6 (t, g), uniform on the set T × G 1 , with upper bounds estimations , respectively, follows from Property A2.Therefore, the absolute boundness of  (25) for real relaxation modulus (36) of the material described by the unimodal Gauss-like distribution, the mean-square identification indices Q N g * N , Equation (24), the mean relative square model approximation index Q Nrel g * N , Equation (38), the sampling points-independent integral indices Q g * N defined by the optimization task (26), and the relative errors ERR (37) of the FMM parameter g * approximation for N measurements independently disturbed by additive, zero mean, normally distributed noises with standard deviation σ = 5 [kPa].

Figure 2 .
Figure 2. Dependence of the parameters of the FMM approximating the "real" relaxation modulu (36): (a)  * ; (b)  * ; (c)  * ; and (d)  * on the number of measurements  for disturbances  = 2, 5, 8 [Pa]; the horizontal purple lines are related to the optimal parameters  * ,  * ,  * and  independent on the sampling instants used in the rheological experiment.

Figure 2 .Figure 3 .
Figure 2. Dependence of the parameters of the FMM approximating the "real" relaxation modulus (36): (a) α * N ; (b) β * N ; (c) G * eN ; and (d) τ * rN on the number of measurements N for disturbances σ = 2, 5, 8 [Pa]; the horizontal purple lines are related to the optimal parameters α * , β * , G * e and τ * r of the experiment data by the model  (,  * ) is confirmed by the values of the relat square model approximation index  ( * ) (38), which for  ≥ 200 and the weak noises does not exceed 0.015%, while for the strongest noises, it does not exceed 0.28 For the noises considered, the values of the model fit indices  ( * ) (24) and  ( (38) and the integral quadratic index ( * ) (26) indicate an excellent fit of the model the experiment data and the fast convergence of  * to  * as  tends to infinity; co pare Tables A1-A3.

Figure 3 .
Figure 3.The indices of the "real" relaxation modulus (36) approximation by the FMM: (a) the mean-square empirical index Q N g * N , Equation (24), (b) the integral quadratic sampling instantsindependent index Q g * N , Equation (26), as the functions of the number of measurements N and noises σ = 2, 5, 8 [Pa]; the horizontal purple lines correspond to the optimal integral index Q(g * ) defined in Equation(27).
For the noises considered, the values of the model fit indices Q N g * N (24) and Q Nrel g * N (38) and the integral quadratic index Q g * N (26) indicate an excellent fit of the model to the experiment data and the fast convergence of g * N to g * as N tends to infinity; compare Tables A1-A3.

Figure 4 .
Figure 4.The measurements G i (red points) of the "real" relaxation modulus (36) and optimal FMM models: sampling points-independent G M (t, g * ) and empirical G M t, g * N for N measurements and normal distribution noises with the standard deviation σ = 8 [Pa]: (a) N = 100; (b) N = 10, 000.

Figure 5 .
Figure 5. Dependence of the mean indices of the "real" relaxation modulus () (36) optimal approximation by the FMM: (a) relative empirical error  (39), (b) integral error  (40) on the number of measurements  and the noises' standard deviations .

Figure 6 .
Figure 6.Dependence of the mean relative error  (41) between the optimal parameters  * and  * of the FMM approximating the "real" relaxation modulus () (36) on the number of measurements  and the noises' standard deviation .

3. 6 .
Material II Consider the material described by the empirical spectrum of relaxation times  introduced by Baumgaertel, Schausberger and Winter [40,41],

Figure 5 .
Figure 5. Dependence of the mean indices of the "real" relaxation modulus G(t) (36) optimal approximation by the FMM: (a) relative empirical error ERRQ Nrel (39), (b) integral error ERRQ (40) on the number of measurements N and the noises' standard deviations σ.

Figure 5 .
Figure 5. Dependence of the mean indices of the "real" relaxation modulus () (36) optimal approximation by the FMM: (a) relative empirical error  (39), (b) integral error  (40) on the number of measurements  and the noises' standard deviations .

Figure 6 .
Figure 6.Dependence of the mean relative error  (41) between the optimal parameters  * and  * of the FMM approximating the "real" relaxation modulus () (36) on the number of measurements  and the noises' standard deviation .

Figure 7 .
Figure 7.The parameters of the FMM approximating the relaxation modulus (43) of mater scribed by the BSW relaxation spectrum (42): (a) * ; (b) * ; (c) * ; and (d) * as the fun of the number of measurements for noises = 3, 6, 8 kPa ; the horizontal purple lines spond to optimal model parameters * , * , * and * being independent on the sampli stants used in the experiment.

Figure 8 .
Figure 8.The indices of the BSW relaxation modulus (42), (43) approximation by the FMM: (a) the mean-square empirical index Q N g * N , Equation (24), (b) the integral quadratic sampling instantsindependent index Q g * N , Equation (26), as the functions of the number of measurements N for the noises σ = 2, 5, 8 [Pa]; the horizontal purple lines correspond to the optimal integral index Q(g * ) defined in Equation(27).

Figure 9 .
Figure 9.The measurements G i (red points) of real relaxation modulus (43) of the material described by the BSW spectrum (42) and the fractional Maxwell optimal models: sampling pointsindependent G M (t, g * ) and empirical G M t, g * N for N measurements and additive random normally distributed noises with standard deviation = 8[Pa] and zero mean value: (a) N = 100; (b) N = 10, 000.

Figure 10 .
Figure 10.Dependence of the mean indices of the "real" BSW relaxation modulus (42), (43) approximation by the FMM: (a) the mean relative empirical error  (39), (b) the mean optimal sampling points-independent integral error  (40) on the number of measurements  and the noises' standard deviation .

Figure 10 .
Figure 10.Dependence of the mean indices of "real" BSW relaxation modulus (42), (43) approximation by the FMM: (a) the mean relative empirical error ERRQ Nrel (39), (b) the mean optimal sampling points-independent integral error ERRQ (40) on the number of measurements N and the noises' standard deviation σ.

Figure 11 .
Figure 11.Dependence of the mean relative error  (41) between the parameters  * and   * of the FMM approximating the "real" BSW relaxation modulus (42), (43) on the number of measurements  and the noises' standard deviation .

Figure 11 .
Figure 11.Dependence of the mean relative error MERR (41) between the parameters g * and g * N of the FMM approximating the "real" BSW relaxation modulus (42), (43) on the number of measurements N and the noises' standard deviation σ.

Appendix A. 1 . 3 .
Uniform on T × G G G 1 Boundness of the FMM Derivative with Respect to β Differentiation of Equation (

Table A2 .
The elements α * N , β * N , G * eN and τ * rN of the FMM parameter vector g * N solving identification task

Table A3 .
(26)elements α * N , β * N , G * eN and τ * rN of the FMM parameter vector g * N solving identification task(25)for real relaxation modulus (36) of the material described by the unimodal Gauss-like distribution, the mean relative square model approximation index Q Nrel g * N , Equation (38), the sampling points-independent integral indices Q g * N defined by the optimization task(26), and the relative errors ERR (37) of the parameter g * approximation for N measurements independently disturbed by additive, zero mean, normally distributed noises with standard deviation σ

Table A4 .
(26)ndix B.2.The Results of the Numerical Studies for Material II For the optimal FMM approximating the relaxation modulus (43) of the material described by the BSW spectrum (42): the elements α * N β * N , G * eN and τ * rN of the vector g * N solving identification task(25), the mean-square identification indices Q N g * N , Equation(24), the mean relative square model approximation index Q Nrel g * N , Equation(38), the sampling points-independent integral indices Q g * N defined by the optimization task(26), and the relative errors ERR (37) of the parameter g * approximation for N relaxation modulus measurements independently disturbed by additive normally distributed noises with standard deviation σ = 3 [kPa].

Table A5 .
(26)the optimal FMM approximating the relaxation modulus(43)of the material described by BSW spectrum (42): the elements α * N , β * N , G * eN and τ * rN of the parameter vector g * N , the meansquare identification indices Q N g * N , Equation(24), the mean relative square model approximation index Q Nrel g * N , Equation(38), the sampling points-independent integral indices Q g * N defined by the optimization task(26), and the relative errors ERR (37) of the parameter g * for N measurements corrupted by the noises with standard deviation σ = 6 [kPa].

Table A6 .
(26)the optimal FMM approximating the relaxation modulus(43)of the material described by the BSW spectrum (42): the elements α * N , β * N , G * eN and τ * rN of the parameter vector g * N , the meansquare identification indices Q N g * N , Equation (24), the mean relative square model approximation index Q Nrel g * N , Equation (38), the sampling points-independent integral indices Q g * N defined by the optimization task(26), and the relative errors ERR (37) of the parameter g * for N measurements corrupted by the noises with standard deviation σ = 8 [kPa].