On Recovery of a Non-Negative Relaxation Spectrum Model from the Stress Relaxation Test Data

The relaxation spectra, from which other material functions used to describe mechanical properties of materials can be uniquely determined, are important for modeling the rheological properties of polymers used in chemistry, food technology, medicine, cosmetics, and many other industries. The spectrum, being not directly accessible by measurement, is recovered from relaxation stress or oscillatory shear data. Only a few models and identification methods take into account the non-negativity of the real spectra. In this paper, the problem of recovery of non-negative definite relaxation spectra from discrete-time noise-corrupted measurements of relaxation modulus obtained in the stress relaxation test is considered. A new hierarchical identification scheme is developed, being applicable both for relaxation time and frequency spectra. Finite-dimensional parametric classes of models are assumed for the relaxation spectra, described by a finite series of power-exponential and square-exponential basis functions. The related models of relaxation modulus are given by compact analytical formula, described by the products of power of time and the modified Bessel functions of the second kind for the time spectrum, and by recurrence formulas based on products of power of time and complementary error functions for frequency spectrum. The basis functions are non-negative. In result, the identification task was reduced to a finite-dimensional linear-quadratic problem with non-negative unknown model parameters. To stabilize the solution, an additional smoothing constraint is introduced. Dual approach was used to solve the stated optimal identification task resulting in the hierarchical two-stage identification scheme. In the first stage, dual problem is solved in two levels and the vector of non-negative model parameters is computed to provide the best fit of the relaxation modulus to experiment data. Next, in second stage, the optimal non-negative spectrum model is determined. A complete scheme of the hierarchical computations is outlined; it can be easily implemented in available computing environments. The model smoothness is analytically studied, and the applicability ranges are numerically examined. The numerical studies have proved that using developed models and algorithm, it is possible to determine non-negative definite unimodal and bimodal relaxation spectra for a wide class of polymers. However, the examples also demonstrated that if the basis functions are non-negative and the model is properly selected for a given type of the real spectrum (unimodal, multimodal), the optimal model determined without non-negativity constraint can be non-negative in the dominant range of its arguments, especially in the wide neighborhood of the spectrum peaks.


Introduction
The viscoelastic relaxation spectrum provides deep insights into the complex behavior of polymers [1][2][3]. The spectrum is not directly measurable and must be recovered from oscillatory shear or relaxation stress data [1,3]. During the last five decades, a number of models and methods have been proposed for the recovery of the relaxation spectrum of a viscoelastic material from oscillatory shear data. The contributions by Baumgaertel and Winter [4], Honerkamp and Weese [5], Malkin [6], Malkin et al. [7], Stadler and Bailly [8], Davis and Goulding [9], Davis et al. [10], and Cho [11] are most frequently cited, as they laid the foundations for several parallel directions of research on the identification of discrete and continuous relaxation spectra based on dynamic modulus data.
Far fewer methods have been proposed for spectrum determination from timemeasurements of the relaxation modulus collected in the stress relaxation test, where time-dependent shear stress is studied for the step increase in the strain. Additionally, some of them only address specific materials. A concise discussion of these works, among which three directions of research can be distinguished, is given in [12]. The three indicated classes of approaches are: (1) differential models and algorithms based on applying the Post-Widder formula [13] of the inverse Laplace transform to designate the relaxation spectrum models proposed in the papers of Alfrey and Doty [14], Ter Haar [15], Bažant and Yunping [16], Goangseup and Bažant [17]; (2) the models derived directly from the known pairs of Laplace transforms proposed by Macey [18], Sips [19,20] and Yamamoto [21] and (3) models based on the expansion of an unknown spectrum into a series of basis functions forming a complete basis in the space of real-valued square-integrable functions developed by Stankiewicz [12,22,23] and Stankiewicz and Gołacki [24]. Some articles are also discussed below.
The relaxation spectra of real materials are non-negative for any relaxation time and any relaxation frequency. However, most of the known models and identification algorithms do not take into account this non-negativity property. Therefore, the resulting spectrum model may take a negative value for some relaxation times or frequencies.
The exceptions are those methods that use the spectrum approximation by non-negative definite simple functions, represented by the Macey [18] exponential-hyperbolic model of the spectrum, the Sips [19,20] model given by difference of two exponential functions, next augmented by Yamamoto [21] to consider long-term modulus. However, resulting spectrum models are positive for all arguments; the scope of their effective applicability is limited due to rather narrow classes of models. The Alfrey and Doty [14] simple differential model, the Ter Harr [15] approximation of the spectrum of relaxation frequencies by the modulus multiplied by the time inverse of the relaxation frequency and other methods using the Post-Widder inversion formula to designate the relaxation spectrum model, as Bažant and Yunping [16] and Goangseup and Bažant [17] two-stage approaches of approximating the stress data by multiple differentiable models of relaxation modulus and next, applying the Post-Widder formula to compute the related spectrum model, guarantee the positive definiteness of the recovered relaxation spectrum whenever the relaxation modulus is a completely monotonic function [25]. Thus, the ranges of their applicability are restricted, also due to the necessity of multiple differentiation of the noise corrupted measurement data. A wider range of applicability has been obtained by Stankiewicz [22] for the non-negative model based on the expanding of an unknown spectrum of relaxation frequencies into a series of basis power-exponential functions. However, article [22] was based on such a definition of the relaxation spectrum, which is not often used in the literature.
Therefore, the goal of the present paper was to formulate and solve the problem of determination of the non-negative definite model of the relaxation spectrum based on discrete-time measurements of the relaxation modulus obtained in the relaxation test.
It was assumed that the approach's proposed and developed identification scheme will be applicable to determine both the relaxation time and frequency spectra. The approximation of the continuous spectrum by finite series of non-negative basis functions was applied. For modeling, the relaxation time spectrum model introduced in [12] was used, while for the spectrum of relaxation frequencies, the basis functions described by the product of power of time and square exponential functions were applied. The components of the relaxation modulus model are given by compact recurrence formulas expressed in

Spectrum of Relaxation
The uniaxial, nonaging, and isothermal stress-strain equation for a linear viscoelastic material can be represented by a Boltzmann superposition integral [3]: where σ(t) and ε(t) denote the stress and strain at the time t and G(t) is the linear relaxation modulus. Modulus G(t) is given by [1,3,12]: or, equivalently, by [1,3]  where H(τ) and H(v) characterize the distributions of relaxation times τ and relaxation frequencies v, respectively. The continuous relaxation spectra H(τ) and H(v), related by H(v) = H 1 v , are generalizations of the discrete Maxwell spectrum [1,3] to a continuous function of the relaxation times τ and frequencies v. Although other definitions of the relaxation spectrum are used in the literature; for example, in [6,22,24,32,33], the definitions introduced by Equations (1) and (2) dominate.
The problem of relaxation spectrum recovery from measurement data, i.e., the problem of solving system of Fredholm integral equations of the first kind (1) or (2), is known to be ill-posed in the Hadamard sense [34], i.e., small changes in measured relaxation modulus can lead to arbitrarily large changes in the determined relaxation spectrum. In remedy, some reduction in the set of admissible solutions can be used. Spectra of relaxation times and frequencies will be modeled by non-negative definite finite series of non-negative basis functions.

Model of Relaxation Time Spectrum
Assume that H(τ) ∈ L 2 (0, ∞), where L 2 (0, ∞) is the space of real-valued squareintegrable functions on the interval (0, ∞). The set of the linearly independent functions e −ατ , τe −ατ , τ 2 e −ατ , . . . form a basis of the space L 2 (0, ∞) [35]; here α is positive timescaling factor. Since the maxima of these basis functions grows rapidly with k, in [12] the scaled basis functions: with the first function are assumed to approximate the relaxation time spectrum H(τ) by the model [12]: where the lower index is the number of model summands. Function h 0 (τ, α) (4) is defined for computational purposes only, since 0 0 = 1 [36], i.e., following [12], the general Formula (3) can be applied in further analysis also for k = 0. Then, according to Equation (1), the respective model of the relaxation modulus is described by: where the basis functions for the spectrum model (6) are given by compact analytical formula specified by the following theorem proved in [12].
Theorem 1 [12]. Let α > 0, k ≥ 0 and t > 0. Then the basis functions ϕ k (t, α) are given by: where K k (x) is the modified Bessel function of the second kind [37] of integer order k.

Model of Relaxation Frequency Spectrum
Assume that the spectrum introduced in Equation (2) is such that H(v) ∈ L 2 (0, ∞). The set of the linearly independent functions e −βv 2 , ve −βv 2 , v 2 e −βv 2 , . . . form a basis of the space L 2 (0, ∞) [38]; here β is a positive time-scaling factor; more precisely, a square of the time-scale factor β expressed in seconds.
Since for any fixed β the maximum: of the function h k (v, β) = v k e −βv 2 , grows or decreases rapidly with k, depending on the value of parameter β, the real relaxation spectrum H(v) can be expanded into a series of normalized basis functions: with the first function as follows where g k are constant model parameters. Since, for the spectra of relaxation times of real materials, the asymptotic property that H(τ) → 0 as τ → ∞ holds, having in mind the relation H(v) = H 1 v , the spectrum of relaxation frequencies tends to zero as the relaxation frequency approaches zero from above, i.e., H(v) → 0 , as v → 0 + . Since = h 0 (0, β) = 1, while = h k (0, β) = 0 for k ≥ 1, the first basis function can be neglected in the series expansion (9). Simultaneously, for practical reasons, it is convenient to replace the infinite summation in (9) with a finite one of K terms, from 1 to K. Therefore, the spectrum H(v) is approximated by a model of the form: where the new basis functions were created as a result of renumbering of = h k (v, β) (8), to unify the presentation of both spectrum models, (5) and (10). According to Equation (2) the respective model of the relaxation modulus G(t) is described by: where The basis functions φ k (t, β) (13) of the model (12) are given by compact recursiveanalytical formulas specified by the following theorem proved in Appendix A.1. Theorem 2. Let β > 0, k ≥ 0 and t ≥ 0. Then the basis functions φ k (t, β) (13) are described by the recursive formula for k ≥ 1, starting with and where the complementary error function er f c(x) is defined by [39]: Thus, the problem of approximating of the continuous spectrum H(v) by finite series H K (v, β) (10) is reduced to problem of the relaxation modulus G(t) approximation by finite linear combination G K (t, β) (12) A few first basis functions h k (v, β) (11) are shown in Figure 1 for two different values of β; the corresponding functions φ k (t, β) (14)- (16) are plotted in Figure 2. It is seen from Figure 1 that the maximum of each scaled basis function h k (v, β) is equal one; however, the relaxation frequency v max yielding to the maximum, for a fixed β, depends on the index k according to the formula v max = k+1 2β , i.e., grows with k. This means that increasing the number of model components K will allow for good modeling of multimodal spectra. However, modeling of such spectra requires a large number of model components, which is confirmed by Example 2 presented below. Reducing the time-scale factor β shifts the spectrum maxima towards larger relaxation frequencies. In turn, from Figure 2, it is seen that the Debye decay monotonicity of basis functions for the relaxation modulus model is in good agreement with the courses of the relaxation modulus obtained in an experiment for real polymers; for example, elastic polyacrylamide hydrogels [28] (Figures 2a,b, 4a, A5, A7, and A8a).
of the relaxation spectrum model ,

Model of Relaxation Frequency Spectrum
Assume that the spectrum introduced in Equation (2) is such that ∈ 0,∞ .
The set of the linearly independent functions , , , … form a basis of the space 0,∞ [38]; here is a positive time-scaling factor; more precisely, a square of the time-scale factor expressed in seconds. Since for any fixed the maximum: ℎ , = of the function ℎ , = , grows or decreases rapidly with , depending on the value of parameter , the real relaxation spectrum can be expanded into a series of normalized basis functions: with the first function ℎ , = , as follows

Positive Definiteness of the Basis Functions
The basis functions of the relaxation frequency spectrum and modulus models are positive definite. Since, for h k (v, β) (11) this property is obvious, the positive definiteness of the functions φ k (t, β) (14)- (16) directly results from their definition (13).

Monotonicity of the Basis Functions
The functions h k (v, β) (11) have a global maximum equal to 1. In view of positive definiteness of basis functions φ k (t, β) (14)- (16), conclusion on their monotonicity results directly from differential property (A1), derived in the Appendix A.1. They are monotonically decreasing for any t ≥ 0.
Despite these properties, in numerical computations, the limited values of φ k (t, β) can be guaranteed only for t ≤ t upp , where t upp depends on the maximal real number accessible in the computing environment. For example, in Matlab the largest finite floatingpoint number in IEEE double precision realmax ∼ = 1.7977·10 308 . Whence, in view of Equation (15), the range of numerical applicability of the model in the time domain is such that e t 2 /4β ≤ realmax, i.e.,

Ranges of Applicability
In models (10) and (12), the parameter β > 0 is a square of the time-scaling factor. The following rule applies: the larger the parameter β, the greater the relaxation times, the lower the relaxation frequencies. The above is illustrated by Figures 1 and 2. Following [12,23], upon the basis functions φ k (t, β) course, the range of applicability is specified as the time t, for which the first K functions φ k (t, β) no longer permanently exceeds, i.e., for any θ > t, ε = 0.5% of its maximum value. Specifically, where Similarly, in [23], the range of applicability specified directly for the relaxation frequencies v was defined based on the variability of the basis functions h The values of t app (β) (23) and v app (β) (24) for different factors β are summarized in Table 1 for K = 5 and K = 12. For K = 6 ÷ 11 the same data are given in Table A1 in Appendix B. Table 1. Ranges of the applicability of the models (10) and (12) for various time-scale parameters for K = 5 and K = 12.

Identification Task
Identification consists in selecting, within the chosen class of models given by Equations (5) and (6) or Equations (10) and (12), of such a model, which ensures the best approximation to the measurement data. To unify the description, we will denote the models G K (t, α) (6) and G K (t, β) (12) of the relaxation modulus, together, as G M (t).
For linear viscoelastic materials, the relaxation modulus is the stress, which is induced in the material when the unit step strain is imposed [3,40]. However, it is impossible to apply a step strain in experiments; loading is never performed infinitely fast [41][42][43]. Therefore, the relaxation modulus is not directly accessible by means of a straightforward measurement method and is usually recovered from the experimental data of the stress relaxation process history collected in non-ideal stress relaxation tests. In such two-phase stress relaxation tests, the strain increases over the loading time interval until a predetermined strain is reached, after which the strain is held constant. Different methods have been proposed during the last few decades for the relaxation modulus determination using the stress data histories from non-ideal relaxation tests [42,[44][45][46][47][48]. The backward recursive method developed by Lee and Knauss [42], the differential rule proposed by Sorvari and Malinen [44], both addressed to the case of constant loading rate, and the general method proposed by Zapas and Phillips [45], where the 'true' relaxation time is delayed of half loading time, are most often cited. For detailed references and an overview, see [41,43,47].
Therefore, suppose a certain identification experiment (stress relaxation test [3,28,40]) resulted in a set of measurements of the relaxation modulus G(t i ) = G(t i ) + z(t i ) at the sampling instants t i ≥ 0, i = 1, . . . , N, where z(t i ) is the measurement noise. It is assumed that the number of measurements N ≥ K. As a measure of the model's accuracy the quadratic index is taken where g K = g 0 · · · g K−1 T is an K-element vector of unknown coefficients of the models (5) and (6) or (10) and (12); · 2 denotes the square norm in the real Euclidean space R N . The matrix Φ N,K is composed of the basis functions ϕ k (t, α) (7) or φ k (t, β) (14)-(16) as follows 26) and G N is the vector of relaxation modulus measurements, i.e., G N = G(t 1 ) · · · G(t N ) T .
For real physical materials, the relaxation spectra H(τ) and H(v) are non-negative for any τ ≥ 0 and v ≥ 0. Thus, the requirement that the respective models H K (τ, α) (5) and H K (v, β) (10) are also non-negative is natural. The basis functions of both classes of models are non-negative. Therefore, if the model parameters are such that g k ≥ 0 for any k = 0, . . . , K − 1, then the models H K (τ, α) and H K (v, β) are non-negative too. The restriction that the model parameters are non-negative is sufficient, but not necessary condition for the non-negativity of the spectrum models. Thus, the optimal identification of non-negative relaxation spectrum models defined by Equations (5) and (6) or Equations (10) and (12) consists in determining the non-negative model parameters minimizing the index Q N (g K , α), i.e., in solving the linear least-squares problem with inequality constraints: where 0 K is K dimensional zero vector. The existence and properties of the solution of the above task depend on the properties of the matrix Φ N,K . Unfortunately, Φ N,K is usually rank-deficient. The linear-quadratic task (27) is ill-conditioned [34] and when the data are noisy, even small changes of the data G N would lead to arbitrarily large artefacts in the optimal model parameters. Therefore, the numerical solution of the finite-dimensional problem (27) is fraught with the same difficulties that the original continuous ill-posed problem of the numerical solution of the Fredholm Equations (1) and (2). The fluctuations of the solution of optimization task (27) may be reduced by introducing an additional direct smoothing constraint g K 2 ≤ κ, where a constant 0 < κ < g N K 2 estimates the level of smoothness assumed for the model parameters g K ; g N K is the normal (minimum Euclidean norm) solution of the original least-squares problem without constraints. As a result, the modified problem of optimal relaxation spectrum identification is obtained: solve minimization task (27) under constraint Dual approach is applied below to solve the optimization task (27), (28).

Results and Discussion
In this section, the optimal identification problem (27) with additional smoothing constraint (28), is solved by applying the dual approach. The existence of the solution of dual maximization task is proved. Next, the idea of parametric optimization [49] is applied to solve the dual task. The necessary and sufficient optimality condition for partial dual tasks is derived in the form of the algebraic polynomial equation. Hierarchical two-stage identification scheme, with the solution of the dual maximization task in two levels, is proposed. Their numerical realization and application of the singular value decomposition technique are discussed. A complete computational algorithm is outlined. The analysis of the smoothness of the relaxation spectra models is presented.

Dual Problem
By introducing Lagrangian multipliers, a vector λ ≥ 0 K and a scalar γ ≥ 0, we can define the Lagrangian for the optimization task (27), (28) where superscript 'T' indicates transpose. The multiplier λ aims at providing a fulfillment of the inequality g K ≥ 0 K . The multiplier γ is the price imposed to satisfy the smoothness constraint (28). The Lagranian is a differentiable function of all arguments. The dual problem takes the form: where the dual function is defined as follows: For an arbitrary κ, λ and γ > 0, the Lagrangian L(g K , λ, γ) is a strictly convex function of g K , which has a unique minimum with respect to g K given by: with symmetric matrix where I K,K is K dimensional identity matrix. Therefore, the dual function defined in Equation (31), by Formulas (29), and (32), after simple algebraic manipulations, can be expressed in compact form as: Before we solve dual problem (30), we will give the algebraic background of the method. The algebraic formula describing L D (λ, γ) will be used to derive the basic result regarding the existence of a solution to the dual problem.

Existence of the Dual Problem Solution
The following proposition, proved in Appendix A.3, is fundamental for the existence of the solution to the optimization task (30). (31), is strictly concave function of both arguments (λ, γ), whenever γ > 0. Now, the existence of the dual problem solution is resolved by the next result proved in Appendix A.4.

Theorem 4.
If the smoothing parameter κ is such that where then there exists a solution λ ,γ of the dual problem (30), such thatγ > 0.

Solution of the Dual Problem
Application of the parametric approach [49] to solve the dual problem (30) results in the scheme: where the function L D (λ, γ) is defined by the following task: From the proof of Theorem 4, it follows that, if the condition (39) is satisfied, then for any fixed λ ≥ 0 K the maximum, with respect to γ, of the strictly concave dual function is positive. Therefore, ∂ ∂γ L D (λ, γ) = 0 is the unique necessary and sufficient condition for γ(λ) optimality in task (42), which, in view of Equations (A19) and (33), immediately yields the following optimality condition for partial dual problem (42).

Theorem 5.
Assume the condition (39) holds. The optimal Lagrange multiplier γ(λ) solves uniquely the optimization task (42) if and only if the following equation is satisfied The above optimality condition means that for any given λ ≥ 0 K and respective optimal γ(λ), the smoothing constraint (28) is satisfied as an equation for the resulting model parameterĝ K (λ, γ(λ)) described by (32). For any fixed λ, Equation (43) can be solved by an arbitrary method of solving nonlinear algebraic equations. By Equations (33), (35) and (36), Equation (43) can be expressed as where, in view of Equation (37), diagonal K × K matrix: (44) can be rewritten as: whence equivalent polynomial equation of unknown variable γ follows: To solve Equation (45) of order 2(r + 1), in general, any numerical method of solving polynomial equations can be applied [51].
In view of Equations (42) and (34) we have: In the Appendix A.5, the following formula is derived: where γ(λ) satisfies Equation (43). By SVD (35) and Equation (36), function L D (λ) (46) takes the form: while the gradient (47), in view of Equation (36), is equivalently given by: For solving the optimization task: numerical methods of constrained nonlinear programming [52] can be applied.

Solution of the Smoothed Identification Problem
If the saddle point of the Lagrangian L(g K , λ, γ) (29) exists, then the dual approach can be successfully applied to solve the optimization task (27) with the smoothing constraint (28), i.e., to solve the stated identification problem. In the case considered, the existence of a saddle point to the Lagrangian follows immediately from Theorem 1, cases (ii) and (iii) in [53] due to the uniqueness of the minimum of L(g K , λ, γ) with respect to g K , given by Equation (32). Thus, the vector: or, equivalently,ĝ whereγ = γ λ , is optimal solution of the optimization task (27), (28), i.e., the vector of the best model parameters. According to Theorem 4, the priceγ > 0 and, by Equation (43), i.e., the smoothness constraint (28) is equally satisfied.

Two-Level Solution od the Dual Problem
To solve the dual problem (30) according to the optimization task (41), i.e., by successive maximization with respect to γ and λ, the following two-level algorithm can be applied.
Upper level of the dual problem solution: Find the multiplierλ ≥ 0 K solving the optimization task (49). The resulting pair λ ,γ ,γ = γ λ , solves the dual problem (30). The numerical computations must be arranged hierarchically, i.e., in each iteration of the maximization procedure at the upper level, the algebraic equation (43) must be solved in the lower-level task. The complete computational procedure for determining the dual problem solution and, next, the optimal model of the relaxation spectrum is given below.

Identification Scheme
The determination of the model of the relaxation spectrum involves the following steps.

1.
Perform the experiment (stress relaxation test [3,28,40]) and record the measurements Choose the number K of model components and, depending on the relaxation spectrum recovery problem considered, the time-scaling factor α for identification of relaxation time spectrum or β for spectrum of relaxation frequencies determination, comparing, for different values of α or β, a few first functions from the sequence Compute the matrix Φ N,K (26) and, next, determine SVD (35). 4.
Compute g N K 2 and choose the constant 0 < κ < g N K 2 .

5.
Determine, in the following two-level computations, the solution λ ,γ of the dual problem (30).

5.1
Choose the initial multiplier λ 0 for the numerical procedure applied to solve optimization task (49).

5.2
Let λ m be the m-th iterate in the numerical procedure solving (49). For λ = λ m solve the Equation (43) according to the chosen numerical procedure and determine γ(λ m ). Polynomial Equation (45) can be solved instead of Equation (43).

5.3
Using γ(λ m ), compute the new multiplier λ m+1 , being the next approximation ofλ, according to the numerical procedure selected to solve the task (49), with the maximized index L D (λ) given by Equation (46) or, equivalently, by Equation (48). If, for λ m+1 , the stopping rule of the chosen numerical procedure is satisfied, i.e., where ε 1 and ε 2 are preselected small positives, takeλ = λ m andγ = γ(λ m ) as the solution to the dual problem (30), and go to Step 6. Otherwise return to Step 5.2 and continue the computations for λ = λ m+1 . 6.
Compute the vector of the optimal model parametersĝ K according to Equation (50) or (51) and, depending on the spectrum recovery problem, determine the optimal model of the relaxation time spectrum given by: or optimal model of the spectrum of relaxation frequencies described by: whereĝ k are elements of the vectorĝ K .
The schematic framework of the above procedure and communication between the two levels solving dual problem and the remaining tasks and relaxation test experiment are shown in Figure 3. or optimal model of the spectrum of relaxation frequencies described by: where are elements of the vector . The schematic framework of the above procedure and communication between the two levels solving dual problem and the remaining tasks and relaxation test experiment are shown in Figure 3.

1.
The SVD (35) of the matrix Φ N,K , of computational complexity O NK 2 [50], must be computed only once in Step 3 and should not be repeated during the two-level computations of Step 5.

2.
The matrix Φ N,K depends on the choice of the basis functions as well as the measurement points t i ; however, it does not depend on the relaxation modulus measurements G(t i ). Thus, when the identification scheme is applied for successive samples of the same material, Step 3 should not be repeated while the same time instants t i are kept and the same model parameters α or β and K are used (selected in Step 2). 3.

5.
In the models proposed the parameters α and β are the time-scaling factors. For the relaxation time model the following rule holds, the lower the parameter α is, the greater the relaxation times are [12]. For the relaxation frequency model, the larger the parameter β, the greater the relaxation times, the lower the relaxation frequencies.
Through the optimal choice of the scaling factors, the best fit of the model to the experimental data can be achieved. In [12], the hierarchical algorithm with the optimal choice of the time-scaling factor α was presented. However, practically in many cases, the selection of the time-scaling factors in Step 2 based on the data concerning model applicability summarized in Tables 1 and A1 for factor β and related tables in [12] for factor α or based on the comparison of a few first functions from the sequences {ϕ k (t, α)} or {φ k (t, β)} for different values of α or β with the experimentally obtained function G(t i ), is quite enough. Similarly, the number K of the models G K (t, α) (6) or G K (t, β) (12) series elements can be initially evaluated. This rough selection strategy of the model parameters selection was applied in [23]. Thus, the choice of K and α must be carried out a posteriori, after the preliminary analysis of the experiment data. 6.
Only the values of λ m and γ(λ m ), not the related parameterĝ K (λ m , γ(λ m )) described by Equation (32), are used in iterations of the numerical procedures solving the dual problem tasks in Step 5. The vectorĝ K of optimal model parameters is computed only in Step 6.

Smoothness Analysis
The smoothing constraint (28) was introduced to stabilize the resulting vectorĝ K (50), for which the equality (52) holds. Since the non-negative basis functions h k (τ, α) (3) and h k (v, β) (11) for any arguments are bounded by one, the following inequalities: hold for the optimal models (53) and (54) with arbitrary time-scale factors, which means that the smoothing of the vector of model parameters imply the boundness of the respective relaxation spectra.
The norms H K (τ, α) 2 and H K (v, β) 2 are also the measures of smoothness of the spectra models, where · 2 means here the square norm in L 2 (0, ∞). Proposition 1 in [Belssel] characterizes H K (τ, α) 2 as the square form of g K with the matrix dependent on α. In the Appendix A.6 the analogous property for the spectrum H K (v, β) is proved.

Proposition 2.
For an arbitrary time-scale factor β and arbitrary vector of model parameters g K for the relaxation spectrum model H K (v, β) (10) we have where Θ is K × K symmetric, positive definite, real matrix of the elements: where k, j = 0, 1, . . . K − 1 and the values of the basis functions φ k (t, β) (14)-(16) at t = 0 are as follows for even indices and for odd indices, where k = 1, 2, . . .. Matrix Θ is a positive definite.
The basis functions φ k (0, 2β) (57), (58) for the time t = 0 do not depend on the time-scale factor, in fact. Propositions 2 and 3 in [12] specify various useful estimates of H K (τ, α) 2 , which can be directly applied to obtain the respective estimates of the norm H K (v, β) 2 and, for the optimal models, results in the next property.
The square roots of the singular values σ 1 (Θ) and σ min (Θ) for K = 5, 6, . . . 12 are summarized in Table 2. Since σ 1 (Θ) grows with K, the greater the number of model summands are, the greater time scaling factor should be, to achieve pre-assumed multiplier σ 1 (Θ)/ 4 2βe in the estimation (59). However, this increase is relatively much weaker than in the case of the model H K (τ, α) (53). Similarly, as in the case of H K (τ, α) (for details, see [12]), estimation (60) is useful only for small K and small time-scale factors. Thus, the smoothness of the vectorĝ K (50) of model parameters guarantees that the fluctuations of the respective spectra of relaxationĤ K (τ, α) andĤ K (v, β) are also bounded. The time-scale factors α and β affect the smoothness of the models.
As in [12], in all examples for numerical experiment N = 5000, sampling instants t i were generated with the constant period in the time interval T = [0, T] seconds and additive measurement noises z(t i ) were selected independently by random choice with uniform distribution on the interval [−0.005, 0.005] Pa. The real spectra and modulus and the basis functions h k (τ, α) (3), h k (v, β) (11) of the spectra models and ϕ k (t, α) (7), φ k (t, β) (14)-(16) of the modulus models were simulated in Matlab R2022a using the special functions besselk and erfc. For the singular value decomposition procedure, svd was applied.
The relaxation time and frequencies models are determined in the class of models defined by H K (τ, α) (5) and H K (v, β) (10). In all examples, the same workflow is applied. First, the optimal models were determined by a two-level regularized least-squares identification scheme proposed in the previous paper [12], i.e., neglecting the non-negativity requirement. This means that, in particular, the optimal time-scale factors α opt and β opt and the optimal regularized model parameter vectors g K were determined resulting in the unconstrained, i.e., determined without non-negativity constraint, models of relaxation time: and relaxation frequencies where g k are elements of the vector g K . Next, for the optimal factors α opt and β opt , the best models with the optimal parameterŝ g K ≥ 0 K are determined with the non-negativity requirement using the scheme proposed above. As a result, the relaxation spectraĤ K (τ, α) (54) andĤ K (v, β) (54) were obtained for time scale factors α = α opt and β = β opt with the non-negative optimal parametersĝ K (50). The smoothing parameter κ was selected several times until a satisfactory accuracy of the fit of the model to the experimental data was obtained. Since some elements of the vectors g K are negative, i.e., the respective components of the models (61) and (62) are negative too, κ smaller than the norm g K 2 are applied.

Example 1
Consider viscoelastic material of relaxation spectrum described by the double-mode Gauss-like distribution considered in [12,27]: inspired by polyethylene data from [27], especially HDPE 1 sample from [27] (Table 1 and Figure 8b), where the parameters are as follows [12]: ϑ 1 = 467 Pa·s, m 1 = 0.0037 s −1 , q 1 = 1.124261 × 10 −6 s −2 , ϑ 2 = 39 Pa·s, m 2 = 0.045 s −1 and q 2 = 1.173 × 10 −3 s −2 . It is shown in [12] that the related real relaxation modulus is Following [12], the time interval T = [0, 1550] seconds is assumed for numerical experiments. In [12], the optimal models H K τ, α opt (61) with the parameter vectors g K and time-scaling factors α opt were determined for K = 3, 4, . . . , 10. Detailed data, including α opt , g K , regularization parameters λ GCV α opt , the square norms g K 2 and H K τ, α opt 2 , as the measures of the solution smoothness, and mean square identification index Q N ( g K )/N were summarized in [12] and (Tables 3 and A3). Only α opt , Q N ( g K )/N and g K 2 are rewritten here in Table 3; the last two to compare with respective data for the constrained non-negative modelĤ K τ, α opt (54). The vectors g K are given in [12] (Table A3), from which it can be seen that some of their elements are negative. For K = 3 one element, for K = 4 two elements, for K = 5, 7, 8, 9 three elements and for K = 6, 10, 11, 12 four elements are negative. Also, the spectra H K τ, α opt (61) are negative for some ranges of the relaxation frequencies, see [12] (Figures 4 and 7a) and also Figure 4, below. The non-negative optimal vectorsĝ K are given in Table A2 in Appendix B. Only two or three elements of these vectors are non-zero, the corresponding elements of the Lagrange multipliers vectors λ are obviously zero. Other numerical data for optimal non-negative modelsĤ K τ, α opt , i.e., square norm g K 2 , identification index Q N (ĝ K )/N and the Lagrange multiplierγ are given in the last columns of Table 3. Figure 4 illustrated the course of the real spectrum H(τ) (63), the unconstrained model H K τ, α opt (61) (blue line) and non-negative model H K τ, α opt (54) (green line) for K = 6, 8, 10, 12. The non-negative modelsĤ K τ, α opt are summarized in Figure 5 for K = 6, 7, . . . , 12. In Figure 6, the models of the relaxation modulus G K t, α opt andĜ K t, α opt corresponding to H K τ, α opt (61) andĤ K τ, α opt (54), respectively, computed according to Equation (6) are plotted for K = 8, 12, where the measurements G(t i ) of the real modulus G(t) (64) are also marked. The optimal models G K t, α opt have been better fitted to the experimental data thanĜ K t, α opt , thus G K t, α opt practically coincide with the measurement points. The deterioration of the identification index for the non-negative models changes from 5.97 times for K = 12 to 65.56 times for K = 5. Table 3. The parameters of the optimal models in Example 1 for the models H K τ, α opt (61) without non-negativity constraint: optimal time-scale factors α opt , the mean quadratic identification indices Q N ( g K )/N (c.f., definition (25)) and the square norms g K 2 for the optimal model parameters g K [12] (Table A3) and for optimal modelsĤ K τ, α opt (54) determined with the non-negativity constraint: the multiplierγ defined by the optimization task (30) and the norms ĝ K 2 (equal to smoothing parameters κ) and identification indices Q N (ĝ K )/N corresponding to non-negative optimal parametersĝ K (50) from Table A2.

Example 2
Consider again the double-mode Gauss-like distribution described by equation (63) . By the formula = ℋ 1⁄ , the respective spectrum of relaxation frequencies is as follows: Figure 6. The measurements G(t i ) of the real relaxation modulus G(t) (64) (red points) from Example 1 and the optimal approximated models determined without non-negativity constraint G K t, α opt and with non-negativity constraintĜ K t, α opt for K model summands: (a) K = 8; (b) K = 12.
The optimal unconstrained models H K v, β opt (62) with the parameter vectors g K and time-scaling factors β opt were determined using the two-level identification scheme proposed in [12] for K = 6, 7, . . . , 21, and β opt , regularization parameters λ GCV β opt , norms g K 2 and mean square identification index Q N ( g K )/N were enclosed in Table 4. The vectors g K are given in Table A3 in Appendix B for selected K, from which it can be seen that the number of negative elements is less than for Example 1. The courses of the unconstrained models H K v, β opt (62) are illustrated by Figure 7, where the real spectrum H(v) (65) and non-negative modelĤ K v, β opt (54) are also given for even K from 6 to 20. For some K the spectra H K v, β opt (62) are negative for some ranges of the relaxation frequencies, see Figure 7b-d for K = 8, 10, 12. For K = 6, the number of model elements is too small to describe the real bimodal spectrum. Only the stronger maximum of the real spectrum is approximated by both models; however, the approximation is more accurate for the non-constrained model, similarly to the approximation of the relaxation modulus measured by identification index. For K = 8, 10, 12, the non-constrained optimal spectrum is negative for some relaxation frequencies, thus applicability of the proposed scheme is necessary to obtain the non-negative model. However, modelĤ K v, β opt is unimodal. For K ≥ 14 the non-constrained spectrum H K v, β opt is non-negative, becomes bimodal and better approximates both the real relaxation modulus and relaxation spectrum than the modelĤ K v, β opt , being still unimodal. Table 4. The parameters of the optimal models in Example 2 for the models H K v, β opt (62) determined without non-negativity constraint: optimal time-scale factors β opt , regularization parameters λ GCV β opt (for details see [12]), the mean quadratic identification indices Q N ( g K )/N (compare (25)) and the norms g K 2 for the optimal model parameters g K given in Table A3 and for optimal modelsĤ K v, β opt (54) determined with the non-negativity constraint: the multiplierγ defined by the optimization task (30) and the norms ĝ K 2 (equal to smoothing parameters κ) and identification indices Q N (ĝ K )/N corresponding to non-negative optimal parametersĝ K (50), see Table A3.

Without Non-Negativity Constraint
With Non-Negativity Constraint    The non-negative optimal vectorsĝ K are given in Table A3. For this model, only two to five elements of these vectors are zero. Other numerical data for optimal non-negative modelsĤ K v, β opt , i.e., norm ĝ K 2 , index Q N (ĝ K )/N and the Lagrange multiplierγ are given in the last columns of Table 4. In Figure 8, the models of the relaxation modulus G K t, β opt andĜ K t, β opt corresponding to H K v, β opt (62) andĤ K v, β opt (54), respectively, computed according to Equation (12) are plotted for K = 12, 20, also the measurements G(t i ) of the real modulus G(t) (64) are given.
For K = 21, vector g K is non-negative, thusĝ K = g K also solves the optimization task (27), (28). The related spectrum modelĤ K v, β opt = H K v, β opt is plotted in Figure 9a, while the relaxation modulus modelĜ K t, β opt corresponding toĤ K v, β opt (54) is given in Figure 9b. The perfect approximation of the relaxation modulus does not match the good approximation of the relaxation spectrum, and the model has also lost its bimodal character. The model already has too many non-zero terms; exactly 21.

Example 3
Now, we consider viscoelastic material of unimodal relaxation spectrum described by distribution: where the parameters are as follows: = 39 • , = 0.045 , and = 1.173 × 10 . The corresponding real relaxation modulus is described by one component of the model of the form (64). In experiment the time interval = [0,500] seconds was applied, which resulted from the inspection of the course of . For = 3,4, … ,8, the optimal time-scaling factors , the related regularization parameters , the mean optimal identification indices / and square norms ‖ ‖ are given in Table 5. The vectors of optimal model parameters are given in Table A4 in Appendix B; the elements of these vectors are both negative and positive. related to model , (62) Next, for time-scale factor = the optimal non-negative models , (54) were determined; the smoothing parameter was selected several times until a satisfactory accuracy of the fit of the model to the experimental data was obtained. The non-negative optimal parameters (50) are given in Table A4, while the multiplies However, for K ≤ 20, as K increases, the models H K v, β opt (62) determined without non-negativity constraint approximates the bimodal spectrum more and more closely, the model determined with the non-negativity constraint does not. ModelsĤ K v, β opt (54), with increasing K, better and better approximate the second, major, maximum, but it does not model the first, minor, maximum even for 20 components.
The corresponding real relaxation modulus G(t) is described by one component of the model of the form (64). In experiment the time interval T = [0, 500] seconds was applied, which resulted from the inspection of the course of G(t). For K = 3, 4, . . . , 8, the optimal time-scaling factors β opt , the related regularization parameters λ GCV , the mean optimal identification indices Q N ( g K )/N and square norms g K 2 are given in Table 5. The vectors of optimal model parameters g K are given in Table A4 in Appendix B; the elements of these vectors are both negative and positive. related to model H K (v, β) (62) Next, for time-scale factor β = β opt the optimal non-negative modelsĤ K (v, β) (54) were determined; the smoothing parameter κ was selected several times until a satisfactory accuracy of the fit of the model to the experimental data was obtained. The non-negative optimal parametersĝ K (50) are given in Table A4, while the multipliesγ defined in Equation (30), norms ĝ 2 , and optimal identification indices Q N (ĝ K )/N are given in Table 5. For K = 3, 4, . . . 12, the real spectrum H(v) (66), optimal models H K v, β opt (62) and H K v, β opt (54) are plotted in Figure 10. The norms ĝ 2 are equal to the smoothing parameters κ, which are assumed smaller than the norms g K 2 , since a quick inspection of the data from Table A4 shows that many elements of the vector g K are negative (even, six by eight for K = 8). In Figure 11 the optimal models of the relaxation modulus G K (t) related to unconstrained H K v, β opt (62) and non-negativeĤ K v, β opt (54) relaxation spectra are plotted for K = 3 and K = 7. Table 5. The parameters of the optimal models in Example 3 for the model H K v, β opt (62) without non-negativity constraint: optimal time-scale factors β opt , regularization parameters λ GCV [12], the mean quadratic identification indices Q N ( g K )/N and the square norms g K 2 for the optimal model parameters g K [12] and for optimal modelĤ K v, β opt (54) determined with the non-negativity constraint: the multiplierγ defined in Equation (30) and the norms ĝ K 2 (equal to smoothing parameters κ) and identification indices Q N (ĝ K )/N corresponding non-negative optimal parameterŝ g K (50).

Without Non-Negativity Constraint
With  Figure 10 shows that the H K v, β opt model, determined without additional nonnegativity constraint, is negative in some range of frequencies for any K. However, the identification index Q N (ĝ K )/N is from 1.08 (for K = 7) to 68.1 (for K = 3) times greater than Q N ( g K )/N obtained for unconstrained model (Table 5), the inspection of Figure 10 shows that modelĤ K v, β opt well approximates the real spectrum, and the quality of this approximation improves with increasing K. For K = 3, . . . , 8, two elements of the vectorĝ K are zero, i.e., for K = 3 only one element of the vectorĝ K is non-zero and in result index Q N (ĝ K )/N is the biggest, see also Figure 11a. Analysis of both the values of identification index Q N (ĝ K )/N and Figure 10 indicates that the best model with non-negativity constraint was obtained for K = 7. For K = 7, both the relaxation modulus models practically coincide with the measurement points and with each other, see Figure 11b. Increasing the number of model components to K = 8 no longer corrects the model.

Discussion
In Example 1, the peaks of the spectrum are more distant than in Example 2. For successive k, the maxima of the basis functions h k (τ, α) (3) are more distant than the maxima of the functions h k (v, β) (11). Thus, the relaxation time model H K (τ, α) (5) was more appropriate for modeling spectrum in Example 1, than model H K (v, β) (10). For the same reason, in Example 1, it was enough to use K = 12 model components, while in Example 2, many more model components were necessary (K = 21) to obtain a satisfactory approximation of the real relaxation modulus and spectrum.
The parameter vectorsĝ K of the models determined with non-negativity constraint have zero elements. Therefore, these models are composed of fewer items than index K would indicate. However, the model of full dimension K must be determined on the identification stage. The proposed approach, effective for the unimodal spectrum, is less effective for the multi-modal spectra, because the additional non-negativity constraint reduces the set of admissible models and, therefore, makes it impossible to achieve such a good fit of the model to the experiment data as for the model determined without this constraint.

Discussion
In Example 1, the peaks of the spectrum are more distant than in Example 2. For successive , the maxima of the basis functions , (3) are more distant than the maxima of the functions ℎ , (11). Thus, the relaxation time model ℋ , was more appropriate for modeling spectrum in Example 1, than model , (10). For the same reason, in Example 1, it was enough to use = 12 model components, while in Example 2, many more model components were necessary ( = 21) to obtain a satisfactory approximation of the real relaxation modulus and spectrum.
The parameter vectors of the models determined with non-negativity constraint have zero elements. Therefore, these models are composed of fewer items than index would indicate. However, the model of full dimension must be determined on the identification stage. The proposed approach, effective for the unimodal spectrum, is less effective for the multi-modal spectra, because the additional non-negativity constraint reduces the set of admissible models and, therefore, makes it impossible to achieve such a good fit of the model to the experiment data as for the model determined without this constraint.
Additionally, the examples showed that a new model of the frequency spectrum can be applied for unimodal and bimodal spectra approximation when the regularized least-squares identification with optimal choice of the time-scale factor is used without additional non-negativity constraint.

Applicability of the Scheme to Discrete Relaxation Spectra Identification
Assume, as above, that the experiment resulted in a set of the measurements ̅ = + at the times ≥ 0, = 1, … , . By (1), for any time we have: Let = ∆ + ∆ , where = 0,1, … , − 1 and ∆ > 0 is the length of integration step. Then, for any = 1, … , , the integral of the right-hand side of Equation (67) can be approximated by: whenever the number of subintervals and the integration step ∆ are such that the integrand is sufficiently small for ≥ − ∆ . Denoting: Figure 11. The measurements G(t i ) (red points) of the real relaxation modulus G(t) from Example 3 and the models of the relaxation modulus corresponding to the optimal spectra models determined without non-negativity constraint H K v, β opt (62) (blue line) and with non-negativity constraint H K v, β opt (54) (green line) for: (a) K = 3; (b) K = 7 summands of the model.

Additionally
, the examples showed that a new model of the frequency spectrum can be applied for unimodal and bimodal spectra approximation when the regularized least-squares identification with optimal choice of the time-scale factor is used without additional non-negativity constraint.

Applicability of the Scheme to Discrete Relaxation Spectra Identification
Assume, as above, that the experiment resulted in a set of the measurements G(t i ) = G(t i ) + z(t i ) at the times t i ≥ 0, i = 1, . . . , N. By (1), for any time t i we have: Let τ k = ∆τ 2 + k∆τ, where k = 0, 1, . . . , K − 1 and ∆τ > 0 is the length of integration step. Then, for any i = 1, . . . , N, the integral of the right-hand side of Equation (67) can be approximated by: whenever the number of subintervals K and the integration step ∆τ are such that the integrand is sufficiently small for τ ≥ K − 1 2 ∆τ. Denoting: compare Equation (26), the set of discretized model equations takes the form: with vector g K of unknown relaxation spectrum at relaxation times τ k and known elements of the matrix Φ N,K , where G M is the vector of the relaxation modulus of model (1) at times t i , defined by analogy to the vector of relaxation modulus measurements G N . Now, the square of the model (70) error is described by identification index Q N (g K ) (25) and the proposed identification scheme can be applied for determining the best nonnegative approximations of the discretized relaxation time spectrum. As a result, the set of pairs τ k ,Ĥ(τ k ) , for k = 0, 1, . . . , K − 1, where the optimalĤ(τ k ) are uniquely given by the optimal model parameterĝ K according to Equation (69). The approximation of the the discrete spectrum becomes more accurate as more rectangles are used in the series (68). By analogous discretization of Equation (2), discrete relaxation frequency spectrum can be determined.
The simple rectangular (midpoint) rule with equally spaced points τ k is applied here; however, other, more sophisticated quadratures can be also used.

Conclusions
In this paper, a new hierarchical identification scheme for recovery of the non-negative continuous relaxation spectra has been derived. The scheme can be applied to identify both relaxation time and frequency spectra using the relaxation test data. Two classes of models are considered; both are based on an expansion of an unknown spectrum into a series of non-negative basis functions. The continuous spectrum of relaxation times was approximated by finite series of power-exponential basis functions, with the components of the relaxation modulus model described by the product of power of time and the modified Bessel function of the second kind. For modeling of the relaxation frequency spectrum, the basis functions described by the product of power of time and square exponential functions were applied. The components of the related relaxation modulus model were proven to be described by compact recurrence formulas expressed in terms of the products of power of time, exponential, and complementary error functions. The quadratic identification index related to the relaxation modulus measurements was used, and an additional smoothing constraint was imposed on the model parameters to ensure the problem was well-posed.
The numerical experiments showed that both considered classes of models can be applied for unimodal and bimodal relaxation spectra modeling with additional non-negativity constraints. The model of the relaxation time spectrum using the modified Bessel functions can be recommended for modeling bimodal spectra with the peaks more distant than the relaxation frequency model. However, the examples showed that in many cases, the non-negative models of the relaxation spectra or models non-negative for almost all arguments can be obtained also using the classical approach, without the additional constraint of the model parameters non-negativity, whenever the basis functions of the relaxation spectrum model are nonnegatively defined. Thus, the following procedure can be recommended for the nonnegative relaxation spectrum determination. First, find the best model of the relaxation spectrum using regularized least-squares identification and check the definiteness of the designated model. If this model is non-negative over a significant range of relaxation times or frequencies, accept it. Otherwise, apply the proposed two-stage hierarchical algorithm and determine the nonnegative relaxation spectrum model. However, the best non-negative model can be obtained by solving the original infinitedimensional task of optimal approximation of the real spectrum in the class of continuous non-negative functions by applying the calculus of variations technique. It will be the subject of further work.  where both the numerator and denominator tends to zero, when t → ∞ . Thus, using the L'Hospital's rule, having in mind (A3), we obtain To prove (20), the product tφ 0 (t, β) is rewritten as , whence, by the L'Hospital's rule we have For k = 1, by (16) and (A10), we immediately have Having in mind (A4), the product tφ 1 (t, β) can be expressed as a fraction , which numerator and denominator tends to zero, when t → ∞ . Thus, the L'Hospital's rule and (A3) yield , whence the next limit directly follows Similarly, to prove (19) for k = 2, the product tφ 2 (t, β), due to (14), (16) and (15), can be expressed as where numerator and denominator tends to zero, when t → ∞ . Thus, using again the L'Hospital's rule and differential formula (A3) we obtain i.e., the dual function grows, as a function of γ, in the near right neighborhood of the zero. To show the above, note that by (A19) and the left inequality in (A20), having in mind that due to (36) and (37) The assumption (39) imply Thus, there exists γ 1 > 0, for example such that for any λ ≥ 0 K , having in mind that also Y K ≥ 0 K , by (A27), the inequality ∂ ∂γ L D (λ, γ) > 0 holds for any 0 < γ < γ 1 , i.e., the dual function increases with growing γ in the near neighborhood of γ = 0. Therefore, the optimal priceγ ≥ γ 1 > 0, which completes the proof.
Appendix B Table A1. Ranges of the applicability of the models (10) and (12) Table A4. Optimal parameters g K andĝ K of the best models H K (v, β) (62) andĤ K (v, β) (54) of the relaxation spectrum from Example 3 determined without and with the non-negativity constraints, respectively, for time scale factors β = β opt given in Table 5; the elements of the vectorsĝ K are expressed in [Pa].
Optimal Model Parameters g K Determined without Non-Negativity Constraint