Two-Level Scheme for Identification of the Relaxation Time Spectrum Using Stress Relaxation Test Data with the Optimal Choice of the Time-Scale Factor

Abstract

The viscoelastic relaxation spectrum is vital for constitutive models and for insight into the mechanical properties of materials, since, from the relaxation spectrum, other material functions used to describe rheological properties can be uniquely determined. The spectrum is not directly accessible via measurement and must be recovered from relaxation stress or oscillatory shear data. This paper deals with the problem of the recovery of the relaxation time spectrum of linear viscoelastic material from discrete-time noise-corrupted measurements of a relaxation modulus obtained in the stress relaxation test. A two-level identification scheme is proposed. In the lower level, the regularized least-square identification combined with generalized cross-validation is used to find the optimal model with an arbitrary time-scale factor. Next, in the upper level, the optimal time-scale factor is determined to provide the best fit of the relaxation modulus to experiment data. The relaxation time spectrum is approximated by a finite series of power–exponential basis functions. The related model of the relaxation modulus is proved to be given by compact analytical formulas as the products of power of time and the modified Bessel functions of the second kind. The proposed approach merges the technique of an expansion of a function into a series of independent basis functions with the least-squares regularized identification and the optimal choice of the time-scale factor. Optimality conditions, approximation error, convergence, noise robustness and model smoothness are studied analytically. Applicability ranges are numerically examined. These studies have proved that using a developed model and algorithm, it is possible to determine the relaxation spectrum model for a wide class of viscoelastic materials. The model is smoothed and noise robust; small model errors are obtained for the optimal time-scale factors. The complete scheme of the hierarchical computations is outlined, which can be easily implemented in available computing environments.

1. Introduction

Numerous mathematical rheological models are used to describe the mechanical properties of viscoelastic materials [1,2,3]. The Maxwell and Kelvin–Voight models are, probably, the best-known rheological models. However, deeper insight into the complex behavior of viscoelastic materials is provided by the relaxation spectrum [2,4,5]. The relaxation spectrum is vital for constitutive models and for insight into the properties of a material, since, from the relaxation spectrum, other material functions used to describe rheological properties can be uniquely determined [4,6,7]. It is commonly used to describe, analyze, compare, and improve the mechanical properties of polymers [4,8,9], concrete [10], asphalt [11], rubber [12], wood [13], glass [14], dough [15], polymeric textile materials [16], geophysics [17] and biological materials [18]. The spectrum is not directly measurable; therefore, it must be recovered from oscillatory shear or stress relaxation data [2,5]. A number of different methods and algorithms have been proposed during the last five decades for the recovery of the relaxation spectrum of a viscoelastic material from oscillatory shear data. These works, using different models, approaches and computational techniques, have included contributions by Baumgaertel and Winter [19], Honerkamp and Weese [20], Malkin [21], Malkin et al. [22], Stadler and Bailly [23], Davis and Goulding [24], Davis et al. [25] and Cho [26]. These studies, but also many others, have created new directions of research on discrete and continuous relaxation spectra identification based on dynamic moduli data, and have been conducted ever since [8,11,27].

However, a classical way of studying viscoelasticity is also the stress relaxation test, where time-dependent shear stress is studied for the step increase in strain [1,2,3,28]. For some materials, for example, highly hydrated biological plants, the stress relaxation test is much easier to perform and is more suitable than dynamic oscillatory shear tests [3,29]. There are only a few papers that deal with spectrum determination from time measurements of the relaxation modulus. Additionally, some of them only address specific materials. The first works came from the turn of the 1940s and 1950s. Alfrey and Doty [30] proposed a simple differential model, being, in fact, the first-order Post–Widder formula [31] for the inverse Laplace transform; for details, see the newer edition [32]. Ter Harr [33] approximated the spectrum of relaxation frequencies by the modulus multiplied by time inverse of the relaxation frequency, which can be thought as the Post–Widder inversion formula of the zero order. After many years, Bažant and Yunping [34] and Goangseup and Bažant [35] introduced a two-stage approach of approximating the stress (or retardation) data via multiple differentiable models of the relaxation modulus (creep compliance) and, next, by applying the Post–Widder formula to designate the related relaxation (retardation) spectrum model. The effectiveness of this approach depends, among other aspects, on the function applied to approximate the relaxation modulus. In [34], a logarithmic–exponential model of the relaxation modulus is proposed, for which the authors state the third-order Post–Widder approximation to be satisfactory.

The relaxation spectrum modeling based on the known pairs of Laplace transforms was initiated by Macey [36], who described the relaxation modulus of viscoelastic ceramic material by the modified Bessel function of the second kind and zero order, which corresponds to the exponential–hyperbolic model of the spectrum. To describe the mechanical properties of polyisobutylene, Sips [37] introduced a simple relaxation spectrum model given by the difference of two exponential functions and a related logarithmic model of the modulus. This model was augmented to consider a long-term modulus by Yamamoto [38] and applied to test the rheological properties of the plant cell wall.

Both the algorithms based on the Post–Widder formula and those using the pairs of Laplace transforms assumed rather narrow classes of models. Thus, the scope of their effective applicability is limited to similar, not much wider classes of ‘real’ relaxation characteristics. A wider range of applicability has been offered by the models based on the expansion of an unknown spectrum into a series of basis functions (or polynomials) forming a complete basis in a function space, for example, in the space of real-valued square-integrable functions. Stankiewicz [39] and Stankiewicz and Gołacki [18] derived algorithms of the optimal regularized identification of relaxation and retardation spectra in the classes of models generated by different special functions, where various rules were applied for the choice of regularization parameters.

However, articles [18,39] were based on such a definition of the relaxation spectrum, according to which, the modulus was directly given by the Laplace integral of the spectrum. This spectrum definition is not often used in the literature. Therefore, computationally efficient algorithms to determine the relaxation spectrum applied to time measurements of the relaxation modulus are still desirable. Recently, for the dominant definition of the relaxation spectrum in the literature, a class of algorithms for relaxation spectrum recovery, which combines the technique of an expansion of a function into a series in an orthonormal basis with the least-squares regularized identification, has been derived by Stankiewicz [40]. Legendre, Laguerre, Chebyshev and Hermite functions were used as the basis functions for the spectrum model. In [40], the problem of determining the spectrum of relaxation frequencies was considered. The complementary problem of determining the spectrum of relaxation times based on the stress relaxation data is both analytically and numerically more difficult than the determination of the spectrum of relaxation frequencies, because the relaxation time occurs in the denominator of the exponential function in the spectrum definitional formula. This very task is the subject of this article. Previous studies [18,39,40] have suggested that by selecting an appropriate time-scale factor, a better fit of the model to the measurement data can be obtained.

The objective of the present paper was to develop a model and an identification algorithm for the determination of the continuous relaxation time spectrum based on discrete-time measurements of the relaxation modulus, which, taking into account the ill-posedness of the original problem of the spectrum recovery and the idea of the optimal selection of the time-scale factor, will provide: (a) good approximation of the relaxation spectrum and modulus also due to the best choice of the time-scale factor; (b) smoothness of the spectrum fluctuations, even for noise-corrupted measurements; (c) noise robustness; (d) applicability to a wide range of viscoelastic materials; (e) ease of the implementation of the model and identification algorithm in available computing packages. The idea of the optimal choice of time-scale factor is used here for the first time in the context of relaxation spectrum identification. The goal of this work was the synthesis of the respective model and identification scheme and the analysis of their properties. A further purpose was the numerical verification of the model and algorithm for double-mode Gauss-like distribution used to describe the viscoelastic properties of various materials, in particular different polymers, wood and glass.

A new hierarchical algorithm of noise robust approximation of the continuous spectrum of relaxation times by finite series of power–exponential basis functions is proposed. The components of the relaxation modulus model are given by compact analytical formulas described by the product of power of time and the modified Bessel function of the second kind. The main properties of the basis functions of relaxation spectrum and modulus models have been studied; positive definiteness, upper bounds, monotonicity and asymptotic properties have been examined. Ranges of applicability for different time-scale factors are determined. Since the problem of relaxation spectrum identification is an ill-posed inverse problem, the regularized least-squares identification technique is applied combined with generalized cross-validation to guarantee the stability of the scheme. The quadratic identification index refers to the measured relaxation modulus. The task of determining the best “regularized” model is solved at the lower level of the identification scheme. The appropriate choice of the time-scale factor on the upper level of the scheme ensures the best fit of the relaxation modulus to experimental data.

The optimality conditions are derived both for the identification problem solved at the lower level and the task of the optimal choice of the time-scale factor at the upper level of the scheme. It is proved that the smoothness of the vector of the optimal model parameters implies smoothness of the fluctuations of the relaxation spectrum model. A direct formula and upper and lower bounds for the square integral norm of the smoothed spectrum model are derived. The accuracy of the spectrum model for noisy measurements of the relaxation modulus is studied, and the linear convergence to the model that we would obtain for the noise-free measurements is proved.

To design a numerical algorithm based on the scheme, the computations should be arranged hierarchically in a two-level scheme, i.e., for each iteration of the minimization procedure at the upper level, the whole numerical procedure must be realized for the lower-level task. The complete computational procedure for determining the best model is described. The singular value decomposition technique is applied to simplify the algebraic computations. The identification scheme can be easily implemented in available computing environments. The numerical studies, especially the strong smoothing of the relaxation spectrum models in the example of a double-mode Gauss-like spectrum, with a simultaneous very good fit to the experimental data of the relaxation modulus models, suggests the need to explore the applicability of the regularized weighted least-squares [41] in the lower-level identification task. This will be the subject of further work.

In Appendix A, the proofs and derivations of some mathematical formulas are given. Some tables have been moved to Appendix B to increase the clarity of the article.

2. Materials and Methods

2.1. Relaxation Time Spectrum

The uniaxial, non-aging and isothermal stress–strain equation for a linear viscoelastic material is represented by a Boltzmann superposition integral [2]:

$$\sigma \left(t\right)={\int }_{-\infty }^{t}G\left(t-u\right)\dot{\epsilon }\left(u\right)du,$$

(1)

where $u$ is the past time variable in the range $-\infty$ to the present time $t$. In Equation (1), variables $\sigma \left(t\right)$ and $\epsilon \left(t\right)$ denote, respectively, the stress and strain at time $t$, and $G\left(t\right)$ is the linear relaxation modulus. Modulus $G\left(t\right)$ is given by [2,5]:

$$G\left(t\right)={\int }_0^{\infty }\frac{H\left(\tau \right)}{\tau }{e}^{-t/\tau }d\tau ,$$

(2)

where $H\left(\tau \right)$ characterizes the distributions of relaxation times $\tau$. The continuous relaxation spectrum $H\left(\tau \right)$ is a generalization of the discrete Maxwell spectrum [2,5] to a continuous function of the relaxation times $\tau$.

The problem of relaxation spectrum identification is the problem of solving a system of Fredholm integral equations of the first kind (2) obtained for discrete-time measurement data. In this paper, time measurements of the relaxation modulus are considered. This problem is ill-posed in the Hadamard sense [42], i.e., small changes in measured relaxation modulus can lead to arbitrarily large changes in the determined relaxation spectrum. As a remedy, some reductions in the set of admissible solutions or the appropriate regularization of the original problem are used. Both the techniques are applied here.

2.2. Models

Assume that $H\left(\tau \right)\in L^2(0,\infty )$, where $L^2(0,\infty )$ is the space of real-valued square-integrable functions on the interval $(0,\infty )$. It is known that the set of the linearly independent functions $\left\{e^{-\alpha \tau },\tau e^{-\alpha \tau },{\tau }^2{e}^{-\alpha \tau },\dots \right\}$ form a basis of the space $L^2(0,\infty )$ [43] (p. 125); here, $\alpha$ is a positive time-scaling factor. Since the maximum

$$\underset{\tau \ge 0}{max} {\overline{h}}_k\left(\tau ,\alpha \right)={\left(\frac{k}{\alpha }\right)}^k{e}^{-k}$$

of the function ${\overline{h}}_k\left(\tau ,\alpha \right)={\tau }^k{e}^{-\alpha \tau }$ grows rapidly with $k$, it is convenient to expand the relaxation spectrum into a series of scaled basis functions

$$h_k\left(\tau ,\alpha \right)={\left(\frac{\alpha \tau }{k}\right)}^k{e}^{-\alpha \tau +k},k=1,2,\dots ,$$

(3)

with the first function

$$h_0\left(\tau ,\alpha \right)=e^{-\alpha \tau },$$

(4)

as follows:

$$H\left(\tau \right)={\sum }_{k=0}^{\infty }{g}_k{h}_k\left(\tau ,\alpha \right),$$

(5)

where $g_k$ are real model coefficients.

It is practical to replace the infinite summation in Equation (5) with a finite one of $K$ first terms, i.e., to approximate the spectrum $H\left(\tau \right)$ by a model of the form

$$H_K\left(\tau ,\alpha \right)={\sum }_{k=0}^{K-1}{g}_k{h}_k\left(\tau ,\alpha \right),$$

(6)

where the lower index is the number of model summands. Then, according to (2), the respective model of the relaxation modulus is described by:

$$G_K\left(t,\alpha \right)={\int }_0^{\infty }\frac{H_K\left(\tau ,\alpha \right)}{\tau }{e}^{-t/\tau }d\tau ={\sum }_{k=0}^{K-1}{g}_k{\varphi }_k\left(t,\alpha \right),$$

(7)

where the functions

$${\varphi }_k\left(t,\alpha \right)={\int }_0^{\infty }\frac{h_k\left(\tau ,\alpha \right)}{\tau }{e}^{-t/\tau }d\tau .$$

(8)

For computational purposes, function $h_0\left(\tau ,\alpha \right)$ is defined by (4). However, since most mathematicians agreed that $0^0=1$ [44], the general Formula (3) can also be applied in further analysis for $k=0$.

The basis functions ${\varphi }_k\left(t,\alpha \right)$ (8) of the modulus model (7) are given by a compact analytical formula specified by the following theorem proved in Appendix A.1.

Theorem 1.

Let $\alpha >0$, $k\ge 0$ and $t>0$. Then, the basis functions ${\varphi }_k\left(t,\alpha \right)$ (8) are given by:

$${\varphi }_k\left(t,\alpha \right)=2e^k{\left(\frac{\sqrt{\alpha t}}{k}\right)}^k{K}_k\left(2\sqrt{\alpha t}\right) ,$$

(9)

where $K_k\left(x\right)$ is the modified Bessel function of the second kind [45,46] of integer order k.

The first Function (9) is as follows

$${\varphi }_0\left(t,\alpha \right)=2K_0\left(2\sqrt{\alpha t}\right) .$$

(10)

The modified Bessel functions of the second kind, and especially $K_0\left(x\right)$ [47], have many applications in science and engineering, for example, in physics to describe the flow of magneto-hydrodynamic (MHD) viscous fluid in a Darcy-type porous medium [48], in engineering to derive a closed analytical form of the model of a axial-flux permanent magnet machine with segmented multipole-Halbach PM array [49] and to describe the per-unit-length internal impedance of two-layer cylindrical conductors [50]. The applications of the Bessel functions in the description of the dynamic response of a mono-pile foundation in homogeneous soil and varied layered soil–rock conditions under horizontal dynamic loads [51], to obtain a fully coupled poroelastic solution for spherical indentation into a half space with an impermeable surface when the indenter is subjected to step displacement loading [52] and to express a distribution of the traveling distance in heterogeneous populations [53] come from material science and ecology.

Five first basis functions $h_k\left(\tau ,\alpha \right)$ (3) are shown in Figure 1 for two different values of the time-scaling factor $\alpha$. Figure 2 shows the related functions ${\varphi }_k\left(t,\alpha \right)$ (9). The basis functions $h_k\left(\tau ,\alpha \right)$ and ${\varphi }_k\left(t,\alpha \right)$ are dimensionless. It is seen from Figure 1 that the maximum of each scaled basis function $h_k\left(\tau ,\alpha \right)$ (3) and (4) is equal one; however, the relaxation time $t^{max}$ corresponding to the maximum, for a given parameter $\alpha$, depends on the index $k$ according to the formula $t^{max}=k/\alpha$; i.e., grows with $k$. This means that increasing the number of model components $K$ will allow for good modeling of multimodal spectra, which is confirmed by the example presented in the final part of the paper. Reducing the time-scale factor $\alpha$ shifts the spectrum maxima towards larger relaxation times. 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 materials; for example: concrete [10] (Figure 13), rubber [12] (Figure 2), elastic polyacrylamide hydrogels [28] (Figures 2a,b, 4a, A5, A7 and A8a), sugar beet [18] (Figure 1) and several foods [3] (Figures 3–39).

2.2.1. Positive Definiteness of the Basis Functions

The basis functions of the relaxation spectrum and modulus models are positive definite. For the functions $h_k\left(\tau ,\alpha \right)$ (3) and (4), these properties are obvious. Since, according to the property A.1 in [54], the Bessel functions of the second kind $K_v\left(x\right)$ are positive for $x>0$ and real $v$, the positive definiteness of the functions ${\varphi }_k\left(t,\alpha \right)$ (9) and ${\varphi }_0\left(t,\alpha \right)$ (10) directly result.

2.2.2. Asymptotic Properties of the Basis Functions

For $x\to 0,$ the following asymptotic formula is found in the literature [45]:

$$K_k\left(x\right)~\frac{1}{2}\Gamma \left(k\right){\left(\frac{x}{2}\right)}^{-k} ,$$

(11)

for modified Bessel functions of the order $k>0$. Thus, for $t\to 0$, Formulas (11) and (9) imply

$${\varphi }_k\left(t,\alpha \right)~\Gamma \left(k\right){\left(\frac{e}{k}\right)}^k ,$$

which means that for $t$ near zero, the values of basis functions ${\varphi }_k\left(t,\alpha \right)$ decrease with increasing index $k$; see Figure 2.

For argument $x\to \infty$, the asymptotic exponential formula holds [54]

$$K_k\left(x\right)~\sqrt{\frac{\pi }{2x}}{e}^{-x} .$$

(12)

From (9) and (12), the next asymptotic formula follows for large $t$ and $k\ge 1$

$${\varphi }_k\left(t,\alpha \right)~\frac{\sqrt{\pi }}{k^k}{\left(\sqrt{\alpha t}\right)}^{k-\frac{1}{4}}{e}^{-2\sqrt{\alpha t}+k} ,$$

(13)

while for $t\to \infty$ and the first basis function, we have

$${\varphi }_0\left(t,\alpha \right)~\sqrt{\frac{\pi }{\sqrt{\alpha t}}}{e}^{-2\sqrt{\alpha t}} .$$

(14)

The multiplication by power function ${\left(\sqrt{\alpha t}\right)}^{k-\frac{1}{4}}$ in (13) causes that the greater $k$ is, the slower the basis function ${\varphi }_k\left(t,\alpha \right)$ decreases, which is seen in Figure 2. The first basis function ${\varphi }_0\left(t,\alpha \right)$ decreases faster than the exponential function, which is also confirmed by the analysis of the course of the basis functions in Figure 2.

By the asymptotic Formulas (13) and (14), the basis functions tend to 0 as $t\to \infty$; i.e., for an arbitrary, $\alpha >0$ and $k=0,1,2,\dots$

$$\underset{t\to \infty }{lim}{\varphi }_k\left(t,\alpha \right)=0.$$

(15)

2.2.3. Upper Bounds for the Basis Functions

In [54], Corollary 3.4, the following inequality is proved:

$$K_0\left(x\right)<\sqrt{\frac{\pi }{2x}} e^{-x},$$

whence

$${\varphi }_0\left(t,\alpha \right)=2K_0\left(2\sqrt{\alpha t}\right)<\sqrt{\frac{\pi }{\sqrt{\alpha t}}} e^{-2\sqrt{\alpha t}} .$$

(16)

From Theorem 3.1 in [54], for $k\ge 1$ we have

$$K_k\left(x\right)<\frac{2^{k-1}\Gamma \left(k\right)}{x^{k+1}}.$$

This inequality applied into (9) gives the following upper bound

$${\varphi }_k\left(t,\alpha \right)<2e^k{\left(\frac{\sqrt{\alpha t}}{k}\right)}^k \frac{2^{k-1}\Gamma \left(k\right)}{{\left(2\sqrt{\alpha t}\right)}^{k+1}}=\frac{e^k}{2k^k} \frac{\Gamma \left(k\right)}{\sqrt{\alpha t}}.$$

(17)

Note that having in mind the positive definiteness of the basis functions ${\varphi }_k\left(t,\alpha \right)$, the limit (15) for the first function ${\varphi }_k\left(t,\alpha \right)$ yields from (16), while, for any $k\ge 1$, can be derived from (17).

2.2.4. Monotonicity of the Basis Functions

As indicated above, the basis functions $h_k\left(\tau ,\alpha \right)$ (3) for $k\ge 1$ have a global maximum equal to 1 for the relaxation time $\tau =k{\alpha }^{-1}$, while the first function $h_0\left(\tau ,\alpha \right)$ (4) is monotonically decreasing.

Since the basis function ${\varphi }_k\left(t,\alpha \right)$ (9) is the product of a monotonically decreasing Bessel function $K_k\left(x\right)$ [54] of increasing argument $2\sqrt{\alpha t}$ and a monotonically increasing power function ${\left(\sqrt{\alpha t}\right)}^k$, its monotonicity is not obvious. By applying the differentiation formula, we obtain [54] (Equation (A.13)):

$$\frac{d}{dx}\left[x^v{K}_v\left(x\right)\right]=-x^v{K}_{v-1}\left(x\right)$$

(18)

which holds for any real-valued order $v$, to (9), gives for integer $k\ge 1$

$$\frac{d}{dt}\left[{\varphi }_k\left(t,\alpha \right)\right]=-4\alpha {\left(\frac{e}{2k}\right)}^k{\left(2\sqrt{\alpha t}\right)}^{k-1}{K}_{k-1}\left(2\sqrt{\alpha t}\right) ,$$

whereas, having in mind the positive definiteness of the Bessel function $K_{k-1}\left(2\sqrt{\alpha t}\right)$ for all $t>0$, we immediately conclude that the function ${\varphi }_k\left(t,\alpha \right)$ is monotonically decreasing; see Figure 2. For the first basis function ${\varphi }_0\left(t,\alpha \right)$ (10), the next differentiation formula is [54] (Equation (A.14)):

$$\frac{d}{dx}\left[K_v\left(x\right)\right]=-\frac{1}{2}\left[K_{v-1}\left(x\right)+K_{v+1}\left(x\right)\right],$$

(19)

which is satisfied for any real $v$, implies

$$\frac{d}{dt}\left[{\varphi }_0\left(t,\alpha \right)\right]=-\sqrt{\frac{\alpha }{t}}\left[K_{-1}\left(2\sqrt{\alpha t}\right)+K_1\left(2\sqrt{\alpha t}\right)\right] ,$$

which shows that ${\varphi }_0\left(t,\alpha \right)$ is also a monotonically decreasing function.

2.2.5. Ranges of Applicability

In the models, the parameter $\alpha >0$ is a time-scaling factor. The following rule holds: the lower the parameter $\alpha$, the greater the relaxation times. The above is illustrated by Figure 1 and Figure 2. Through the optimal choice of the scaling factor, the best fit of the model to the experimental data is achieved, which will be the subject of study in Section 3.2.

Following [40], on the basis of the relaxation modulus course, the range of applicability is specified as the time $t$, for which the first $K$ basis functions ${\varphi }_k\left(t,\alpha \right)$ no longer permanently exceed; i.e., for any $\theta >t$, $\epsilon =0.5\%$ of its maximum value. Specifically,

$$t_{app}\left(\alpha \right)=\underset{0\le k\le K-1}{max}\underset{t>0}{min} \left\{t:\left|{\varphi }_k\left(\theta ,\alpha \right)\right|\le 0.005·{\varphi }_{kmax}\left(\alpha \right)\mathrm{for}\mathrm{any}\theta \ge t\right\},$$

(20)

where

$${\varphi }_{kmax}\left(\alpha \right)=\underset{t\ge 0}{max} \left|{\varphi }_k\left(t,\alpha \right)\right|.$$

Similarly, in [40], the range of applicability specified directly for the relaxation times $\tau$ was defined on the basis of the variability in the basis functions $h_k\left(\tau \right)$; i.e.,

$${\tau }_{app}\left(\alpha \right)=\underset{0\le k\le K-1}{max}\underset{\tau >0}{min} \left\{\tau :\left|h_k\left(\vartheta ,\alpha \right)\right|\le 0.005·h_{kmax}\left(\alpha \right)\mathrm{for}\mathrm{any}\vartheta \ge \tau \right\},$$

(21)

with $h_{kmax}\left(\alpha \right)$ defined by

$$h_{kmax}\left(\alpha \right)=\underset{\tau \ge 0}{max} \left|h_k\left(\tau ,\alpha \right)\right|.$$

The times $t_{app}\left(\alpha \right)$ (20) and ${\tau }_{app}\left(\alpha \right)$ (21) for different values of $\alpha$ are summarized in

Table 1. Ranges of the applicability of the model for various time-scale parameters for $K=5$ and $K=12$.

Time-Scale Factor α [s]	$\mathit{K}=5$	$\mathit{K}=5$	$\mathit{K}=12$	$\mathit{K}=12$
	Range 1 of Relaxation Times τapp(α) [s]	Range 1 of Times tapp(α) [s]	Range 1 of Relaxation Times τapp(α) [s]	Range 1 of Times tapp(α) [s]
0.0001	144,305.22	282,360.6	255,824.35	662,077.1
0.001	14,430.52	28,236.06	25,582.435	66,207.71
0.01	1443.05	2823.61	2558.243	6620.77
0.1	144.305	282.36	255.824	662.08
1	14.43	28.236	25.5824	66.208
10	1.4431	2.824	2.5582	6.621
100	0.1443	0.282	0.2558	0.662

¹ The upper bounds $t_{app}$ (20) and ${\tau }_{app}\left(\alpha \right)$ (21) of the applicability intervals $\left[0,t_{app}\left(\alpha \right)\right]$ and $\left[0,{\tau }_{app}\left(\alpha \right)\right]$ are given.

for $K=5$ and $K=12$. The same data for $K=6,\dots ,11$ are given in

Table A1. Ranges of the applicability of the model for various time-scale parameters for $K=6,\dots 11$.

Time-Scale Factor α [s]	Range 1 of Relaxation Times τapp(α) [s]	Range 1 of Times tapp(α) [s]	Range 1 of Relaxation Times τapp(α) [s]	Range 1 of Times tapp(α) [s]
	$\mathit{K}=6$	$\mathit{K}=6$	$\mathit{K}=7$	$\mathit{K}=7$
0.0001	161,655.35	337,647.7	178,346.044	392,372
0.001	16,165.53	33,764.77	17,834.606	39,238.5
0.01	1616.55	3376.48	1783.462	3923.85
0.1	161.65	337.65	178.348	392.5
1	16.165	33.765	17.836	39.26
10	1.616	3.376	1.784	3.95
100	0.1616	0.337	0.1783	0.392
$\mathit{\alpha }$ $\left[\mathrm{s}\right]$	$K=8$	$K=8$	$K=9$	$K=9$
0.0001	194,529.30	446,719.4	210,306.21	500,804.6
0.001	19,452.93	44,671.94	21,030.621	50,080.46
0.01	1945.29	4467.19	2103.062	5008.046
0.1	194.53	446.72	210.306	500.81
1	19.453	44.672	21.0306	50.08
10	1.945	4.467	2.1031	5.008
100	0.1945	0.447	0.2103	0.501
$\mathit{\alpha }$ $\left[\mathrm{s}\right]$	$K=10$	$K=10$	$K=11$	$K=11$
0.0001	225,748.07	554,697	240,907.43	608,443.6
0.001	22,574.807	55,469.7	24,090.743	60,844.38
0.01	2257.481	5546.97	2409.074	6084.45
0.1	225.748	554.69	240.907	608.46
1	22.5748	55.469	24.0907	60.86
10	2.2575	5.547	2.4091	6.1
100	0.2257	0.555	0.2409	0.608

¹ The upper bounds $t_{app}\left(\alpha \right)$ (20) and ${\tau }_{app}\left(\alpha \right)$ (21) of the applicability intervals $\left[0,t_{app}\left(\alpha \right)\right]$ and $\left[0,{\tau }_{app}\left(\alpha \right)\right]$ are given.

in Appendix B.

A review of the data from these tables shows that for any fixed $\alpha$, both ${\tau }_{app}\left(\alpha \right)$ and $t_{app}\left(\alpha \right)$ grow almost linearly with the number of model summands $K$.

2.3. Least-Squares Regularized Identification

Identification consists of the selection, within the chosen class of models given by (6) and (7), of such a model that ensures the best fit to the measurement results. Suppose a certain identification experiment (stress relaxation test [2,3]) resulted in a set of measurements of the relaxation modulus $\left\{\overline{G}\left(t_i\right)=G\left(t_i\right)+z\left(t_i\right)\right\}$ at the sampling instants $t_i\ge 0$, $i=1,\dots ,N$, where $z\left(t_i\right)$ is the measurement noise. It is assumed that the number of measurements $N\ge K$. As a measure of the model (7) accuracy, the quadratic index is taken

$$Q_N\left({\mathit{g}}_K,\alpha \right)={\sum }_{i=1}^N{\left[\overline{G}\left(t_i\right)-G_K\left(t_i,\alpha \right)\right]}^2,$$

(22)

where ${\mathit{g}}_K={\left[\begin{array}{ccc}{g}_0& \cdots & g_{K-1}\end{array}\right]}^T$ is an $K$-element vector of unknown coefficients of the models (6) and (7). The identification index (22) is rewritten in the compact form as

$$Q_N\left({\mathit{g}}_K,\alpha \right)={\Vert {\overline{\mathit{G}}}_N-{\mathit{\Phi }}_{N,K}\left(\alpha \right){\mathit{g}}_K\Vert }_2^2,$$

(23)

where

$${\mathit{\Phi }}_{N,K}\left(\alpha \right)=\left[\begin{array}{ccc}{\varphi }_0\left(t_1,\alpha \right)& \cdots & {\varphi }_{K-1}\left(t_1,\alpha \right)\\ ⋮& \ddots & ⋮\\ {\varphi }_0\left(t_N,\alpha \right)& \cdots & {\varphi }_{K-1}\left(t_N,\alpha \right)\end{array}\right],{\overline{\mathit{G}}}_N=\left[\begin{array}{c}\overline{G}\left(t_1\right)\\ ⋮\\ \overline{G}\left(t_N\right)\end{array}\right]$$

(24)

and ${\Vert ·\Vert }_2$ denotes the square norm in the real Euclidean space ${ℛ}^N$. Thus, the optimal identification of the relaxation spectrum in the class of models defined by (6) and (7) consists of determining the model parameter ${\mathit{g}}_K$ through solving the following optimization task:

$$\underset{{\mathit{g}}_K\in {ℛ}^K}{min} {\Vert {\overline{\mathit{G}}}_N-{\mathit{\Phi }}_{N,K}\left(\alpha \right) {\mathit{g}}_K\Vert }_2^2.$$

(25)

The matrix ${\mathit{\Phi }}_{N,K}\left(\alpha \right)$ is usually ill-conditioned. In consequence, the problem (25) is ill posed in the sense of Hadamard [42]. The parameter ${\mathit{g}}_K$ minimizing identification index (23) is not unique, and even the normal (with the lowest Euclidean norm) solution of (25) is a non-continuous and unbounded function of the measurement vector ${\overline{\mathit{G}}}_N$. Therefore, when the data are noisy, even small changes in ${\overline{\mathit{G}}}_N$ would lead to arbitrarily large artefacts in optimal model parameters ${\mathit{g}}_K$. To deal with the ill-posedness, Tikhonov regularization [55], replacing the original ill-posed problem (25) by a nearby problem with a modified square functional of the form:

$$\underset{{\mathit{g}}_K\in {ℛ}^K}{min}{\Vert {\overline{\mathit{G}}}_N-{\mathit{\Phi }}_{N,K}\left(\alpha \right) {\mathit{g}}_K\Vert }_2^2+\lambda {\Vert {\mathit{g}}_K\Vert }_2^2,$$

(26)

where $\lambda >0$ is a regularization parameter, can be used. The regularized task (26) is well-posed; that is, the solution always exists, is unique and continuously depends on both the matrix ${\mathit{\Phi }}_{N,K}$ and on the measurement data ${\overline{\mathit{G}}}_N$. The parameter vector solving (26) is given by:

$${\overline{\mathit{g}}}_K^{\lambda }\left(\alpha \right)={\left({\mathit{\Phi }}_{N,K}^T\left(\alpha \right){\mathit{\Phi }}_{N,K}\left(\alpha \right)+\lambda {\mathbb{I}}_{K,K}\right)}^{-1}{\mathit{\Phi }}_{N,K}^T\left(\alpha \right){\overline{\mathit{G}}}_N,$$

(27)

where ${\mathbb{I}}_{K,K}$ is the $K$-dimensional identity matrix.

Following [40], the generalized cross-validation GCV [42,56] is applied, according to which, the best regularization parameter is [42,56]:

$${\lambda }_{GCV}\left(\alpha \right)=min\left\{\lambda : \lambda =arg \underset{\lambda \ge 0}{min} V_{GCV}\left(\lambda ,\alpha \right)\right\},$$

(28)

where the GCV functional is defined by [56]

$$V_{GCV}\left(\lambda ,\alpha \right)=\frac{{\Vert \mathit{\Xi }\left(\lambda ,\alpha \right){\overline{\mathit{G}}}_N\Vert }_2^2}{tr{\left[\mathit{\Xi }\left(\lambda ,\alpha \right)\right]}^2},$$

(29)

with the matrix

$$\mathit{\Xi }\left(\lambda ,\alpha \right)={\mathbb{I}}_{N,N}-{\mathit{\Phi }}_{N,K}\left(\alpha \right){\left({\mathit{\Phi }}_{N,K}^T\left(\alpha \right){\mathit{\Phi }}_{N,K}\left(\alpha \right)+\lambda {\mathbb{I}}_{K,K}\right)}^{-1}{\mathit{\Phi }}_{N,K}^T\left(\alpha \right).$$

(30)

Here, $\mathit{\Xi }\left(\lambda ,\alpha \right){\overline{\mathit{G}}}_N$ is the residual vector for the regularized solution (27); $tr\left[\mathit{\Xi }\left(\lambda ,\alpha \right)\right]$ denotes the trace of the symmetric matrix $\mathit{\Xi }\left(\lambda ,\alpha \right)$. The optimization problem in (28) has a unique solution, and the resulting parameter ${\overline{\mathit{g}}}_K^{{\lambda }_{GCV}\left(\alpha \right)}\left(\alpha \right)$ differs the least from the normal solution of the problem (26) that we would obtain for the ideal (not noise-corrupted) measurements of the relaxation modulus [56].

3. Results and Discussion

In this section, the necessary optimality condition for the regularized identification task with GCV’s choice of regularization parameter is given in the form of an algebraic nonlinear equation. Next, the problem of the choice of the optimal scale factor is stated, and the respective necessary optimality condition is derived. A hierarchical two-level identification scheme with the optimal choice of the scale factor is proposed. The numerical realization and the application of a singular value decomposition technique are discussed. A complete computational procedure is outlined. The analysis of the smoothing of the model and model accuracy for noisy measurements of the relaxation modulus is presented.

3.1. Necessary Optimality Condition

The GCV functional (29) is a differentiable function of the regularization parameter $\lambda$. The necessary optimality condition for the minimization task solved in (28) is derived in Appendix A.3.

Theorem 2.

The optimal regularization parameter ${\lambda }_{GCV}\left(\alpha \right)$ solves the following equation:

$$\begin{gathered} {\overline{\mathit{G}}}_N^T{\left(\mathit{\Psi }\left(\alpha \right)+\lambda {\mathbb{I}}_{N,N}\right)}^{-3}{\overline{\mathit{G}}}_{N}tr\left[{\left(\mathit{\Psi }\left(\alpha \right)+\lambda {\mathbb{I}}_{N,N}\right)}^{-1}\right]={\overline{\mathit{G}}}_N^T(\mathit{\Psi }\left(\alpha \right)+ \\ \lambda {\mathbb{I}}_{N,N}){}^{-2}{\overline{\mathit{G}}}_N tr\left[{\left(\mathit{\Psi }\left(\alpha \right)+\lambda {\mathbb{I}}_{N,N}\right)}^{-2}\right]. \end{gathered}$$

(31)

where

$$\mathit{\Psi }\left(\alpha \right)={\mathit{\Phi }}_{N,K}\left(\alpha \right){\mathit{\Phi }}_{N,K}^T\left(\alpha \right).$$

(32)

In Appendix A.4, the following formula:

$$\begin{gathered} \frac{\partial V_{GCV}\left(\lambda ,\alpha \right)}{\partial \lambda }= \\ \frac{-2{\overline{\mathit{G}}}_N^T{\left(\mathit{\Psi }\left(\alpha \right)+\lambda {\mathbb{I}}_{N,N}\right)}^{-3} {\overline{\mathit{G}}}_{N}tr\left[{\left(\mathit{\Psi }\left(\alpha \right)+\lambda {\mathbb{I}}_{N,N}\right)}^{-1}\right]+2{\overline{\mathit{G}}}_N^T{\left(\mathit{\Psi }\left(\alpha \right)+\lambda {\mathbb{I}}_{N,N}\right)}^{-2}{\overline{\mathit{G}}}_{N}tr\left[{\left(\mathit{\Psi }\left(\alpha \right)+\lambda {\mathbb{I}}_{N,N}\right)}^{-2}\right]}{{\left[tr\left[{\left(\mathit{\Psi }\left(\alpha \right)+\lambda {\mathbb{I}}_{N,N}\right)}^{-1}\right]\right]}^3} \end{gathered}$$

(33)

is developed, describing the derivative of the GCV functional (29) minimized in (28) directly as a function of matrices ${\overline{\mathit{G}}}_N$ and $\mathit{\Psi }\left(\alpha \right)$.

3.2. Choice of the Time-Scale Factor

The optimal regularization parameter ${\lambda }_{GCV}\left(\alpha \right)$ and matrix ${\mathit{\Phi }}_{N,K}\left(\alpha \right)$ depend on the time-scale factor. Therefore, the optimal model parameter

$${\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left(\alpha \right)={{\overline{\mathit{g}}}_K^{\lambda }\left(\alpha \right)|}_{\lambda ={\lambda }_{GCV}\left(\alpha \right)}$$

(34)

and the optimal identification index

$$Q_{Nopt}\left(\alpha \right)=Q_N\left({\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left(\alpha \right),\alpha \right)={\Vert {\overline{\mathit{G}}}_N-{\mathit{\Phi }}_{N,K}\left(\alpha \right){\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left(\alpha \right)\Vert }_2^2$$

(35)

also depend on $\alpha$. By the optimal choice of the scaling factor, the best fit of the model to the experimental data can be achieved. By (27), (30) and (34), and having (A9) in mind, the index $Q_{Nopt}\left(\alpha \right)$ (35) is expressed by the following formula:

$$Q_{Nopt}\left(\alpha \right)={\lambda }_{GCV}^2\left(\alpha \right){\overline{\mathit{G}}}_N^T{\left(\mathit{\Psi }\left(\alpha \right)+{\lambda }_{GCV}\left(\alpha \right){\mathbb{I}}_{N,N}\right)}^{-2}{\overline{\mathit{G}}}_N.$$

(36)

The smaller the index $Q_{Nopt}\left(\alpha \right)$, the smaller the model error results. Thus, the problem of the choice of the best time-scale factor ${\alpha }_{opt}$ takes the form

$$\underset{\alpha >0}{min}{Q}_{Nopt}\left(\alpha \right)=Q_{Nopt}\left({\alpha }_{opt}\right).$$

(37)

Before we state the necessary optimality condition of the task (37), we prove the following result concerning the basis functions ${\varphi }_k\left(t,\alpha \right)$ as the functions of the time-scale factor. The proof is given in Appendix A.5.

Theorem 3.

Let $\alpha >0$ and $t>0$. Then, for $k\ge 1$, derivatives of the basis functions ${\varphi }_k\left(t,\alpha \right)$ (9) with respect to the parameter $\alpha$ are given by:

$$\frac{d}{d\alpha }\left[{\varphi }_k\left(t,\alpha \right)\right]=-e\left(\frac{t}{k}\right){\left(\frac{k-1}{k}\right)}^{k-1}{\varphi }_{k-1}\left(t,\alpha \right).$$

(38)

For $k=0,$ the following formula holds:

$$\frac{d}{d\alpha }\left[{\varphi }_0\left(t,\alpha \right)\right]=-\frac{1}{e\alpha }{\varphi }_1\left(t,\alpha \right) .$$

(39)

Thus, positive definite functions ${\varphi }_k\left(t,\alpha \right)$ are monotonically decreasing functions of $\alpha$ for any fixed $t>0$.

The following necessary optimality condition is derived in Appendix A.6:

Theorem 4.

The optimal time-scale factor ${\alpha }_{opt}$ satisfies the following equation:

$$\begin{gathered} {\lambda }_{GCV}\left(\alpha \right){\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-1}{\mathit{\Omega }}_{N,K}\left(\alpha \right)\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}{\overline{\mathit{G}}}_N=\frac{d{\lambda }_{GCV}\left(\alpha \right)}{d\alpha }{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}{\overline{\mathit{G}}}_N- \\ {\lambda }_{GCV}\left(\alpha \right)\frac{d{\lambda }_{GCV}\left(\alpha \right)}{d\alpha }{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-3}{\overline{\mathit{G}}}_N, \end{gathered}$$

(40)

where derivative $\frac{d{\lambda }_{GCV}\left(\alpha \right)}{d\alpha }$ is given by the equation

$$\begin{gathered} \frac{d{\lambda }_{GCV}\left(\alpha \right)}{d\alpha }\{3{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-4}{\overline{\mathit{G}}}_{N}tr\left[\mathit{{\rm Y}}{\left(\alpha \right)}^{-1}\right]-{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-3}{\overline{\mathit{G}}}_{N}tr\left[\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}\right] \\ -2{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}{\overline{\mathit{G}}}_N tr\left[\mathit{{\rm Y}}{\left(\alpha \right)}^{-3}\right]\}= \\ 2{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}{\mathit{\Omega }}_{N,K}\left(\alpha \right)\mathit{{\rm Y}}{\left(\alpha \right)}^{-1}{\overline{\mathit{G}}}_N tr\left[\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}\right]+ \\ 2{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}{\overline{\mathit{G}}}_N tr\left[{\mathit{\Omega }}_{N,K}\left(\alpha \right)\mathit{{\rm Y}}{\left(\alpha \right)}^{-3}\right]- \\ 2{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-3}{\mathit{\Omega }}_{N,K}\left(\alpha \right)\mathit{{\rm Y}}{\left(\alpha \right)}^{-1}{\overline{\mathit{G}}}_{N}tr\left[\mathit{{\rm Y}}{\left(\alpha \right)}^{-1}\right]- \\ {\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}{\mathit{\Omega }}_{N,K}\left(\alpha \right)\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}{\overline{\mathit{G}}}_{N}tr\left[\mathit{{\rm Y}}{\left(\alpha \right)}^{-1}\right]-{\overline{\mathit{G}}}_N^T\mathit{{\rm Y}}{\left(\alpha \right)}^{-3}{\overline{\mathit{G}}}_{N}tr\left[{\mathit{\Omega }}_{N,K}\left(\alpha \right)\mathit{{\rm Y}}{\left(\alpha \right)}^{-2}\right], \end{gathered}$$

(41)

with symmetric matrices

$$\mathit{{\rm Y}}\left(\alpha \right)=\mathit{\Psi }\left(\alpha \right)+{\lambda }_{GCV}\left(\alpha \right){\mathbb{I}}_{N,N},$$

(42)

$${\mathit{\Omega }}_{N,K}\left(\alpha \right)=\frac{d}{d\alpha }\mathit{\Psi }\left(\alpha \right)={\mathit{\Theta }}_{N,K}\left(\alpha \right){\mathit{\Phi }}_{N,K}^T\left(\alpha \right)+{\mathit{\Phi }}_{N,K}\left(\alpha \right){\mathit{\Theta }}_{N,K}^T\left(\alpha \right).$$

(43)

where $\mathit{\Psi }\left(\alpha \right)$ is defined by (32) and

$$\begin{gathered} {\mathit{\Theta }}_{N,K}\left(\alpha \right)=\frac{d}{d\alpha }{\mathit{\Phi }}_{N,K}\left(\alpha \right)= \\ -\left[\begin{array}{cccc}\frac{1}{e\alpha }{\varphi }_1\left(t_1,\alpha \right)& e{t}_1{\varphi }_0\left(t_1,\alpha \right)& \cdots & e\left(\frac{t_1}{K-1}\right){\left(\frac{K-2}{K-1}\right)}^{K-2}{\varphi }_{K-2}\left(t_1,\alpha \right)\\ ⋮& ⋮& \ddots & ⋮\\ \frac{1}{e\alpha }{\varphi }_1\left(t_N,\alpha \right)& e{t}_N{\varphi }_0\left(t_N,\alpha \right)& \cdots & e\left(\frac{t_N}{K-1}\right){\left(\frac{K-2}{K-1}\right)}^{K-2}{\varphi }_{K-2}\left(t_N,\alpha \right)\end{array}\right]. \end{gathered}$$

(44)

3.3. Two-Level Identification Scheme

To find the optimal time-scale factor ${\alpha }_{opt}$ and the optimal model of the relaxation time spectrum, the following two-level scheme is applied:

3.3.1. Lower Level

Given time-scale factor $\alpha >0$, find regularization parameter ${\lambda }_{GCV}\left(\alpha \right)$ solving the GCV minimization task (28).

3.3.2. Upper Level

Find the time-scale factor ${\alpha }_{opt}>0$ minimizing identification index $Q_{Nopt}\left(\alpha \right)$ (36), i.e., solving optimization task (37). Take

$${\widehat{\mathit{g}}}_K={\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }_{opt}\right),$$

(45)

where ${\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }_{opt}\right)$ is defined by (34) as a parameter of the best model of the relaxation time spectrum.

Having the optimal model parameter ${\widehat{\mathit{g}}}_K$, the optimal model of the relaxation spectrum is determined according to the formula resulting directly from (6):

$$H_{Kopt}\left(\tau \right)=H_K\left(\tau ,{\alpha }_{opt}\right)={\sum }_{k=0}^{K-1}{\widehat{g}}_k{h}_k\left(\tau ,{\alpha }_{opt}\right),$$

(46)

where ${\widehat{g}}_k$ are elements of the vector ${\widehat{\mathit{g}}}_K$.

To design numerical realization of the scheme, we need:

• A numerical procedure for solving the lower-level GCV minimization task (28);
• An iterative scheme for solving the upper-level problem (37) of choosing the best time-scale factor.

For any given parameter $\alpha$, the GCV function $V_{GCV}\left(\lambda ,\alpha \right)$ (29) is differentiable with respect to $\lambda$. Partial derivative $\frac{\partial V_{GCV}\left(\lambda ,\alpha \right)}{\partial \lambda }$ is given by (33) as a function of the experiment data ${\overline{\mathit{G}}}_N$ and the matrix $\mathit{\Psi }\left(\alpha \right)$ (32), which depends on the time instants $t_i$ used in the relaxation experiment and on the time-scaling factor. An arbitrary gradient optimization method can be implemented to solve the GCV minimization task (28). Additionally, an arbitrary gradient method can be used to solve optimization problem (37); the derivative of the index $Q_{Nopt}\left(\alpha \right)$ (36) is described by Formula (A24) derived in Appendix A.6. However, the optimal parameter ${\alpha }_{opt}$ can be also found by solving the necessary optimality condition from Theorem 4, i.e., the two scalar algebraic Equations (40) and (41), in fact.

3.4. Algebraic Background of the Identification Scheme

Formulas (27)–(30), fundamental for lower-level optimization task (28), are elegant but generally unsuitable for computational purposes. Following [40], the singular value decomposition (SVD) technique [57] will be applied. Let SVD of the $N\times K$ dimensional matrix ${\mathit{\Phi }}_{N,K}\left(\alpha \right)$ take the form [57]:

$${\mathit{\Phi }}_{N,K}\left(\alpha \right)=\mathit{U}\left(\alpha \right)\mathit{\Sigma }\left(\alpha \right) \mathit{V}{\left(\alpha \right)}^T,$$

(47)

where $\mathit{\Sigma }\left(\alpha \right)=diag\left({\sigma }_1\left(\alpha \right),\dots ,{\sigma }_{r\left(\alpha \right)}\left(\alpha \right),0,\dots ,0\right)ϵ{ℛ}^{N,K}$ is a diagonal matrix containing the non-zero singular values ${\sigma }_1\left(\alpha \right),\dots ,{\sigma }_{r\left(\alpha \right)}\left(\alpha \right)$ of the matrix ${\mathit{\Phi }}_{N,K}\left(\alpha \right)$, matrices $\mathit{V}\left(\alpha \right)\in {ℛ}^{K,K}$ and $\mathit{U}\left(\alpha \right)\in {ℛ}^{N,N}$ are orthogonal and $r\left(\alpha \right)=rank\left[{\mathit{\Phi }}_{N,K}\left(\alpha \right)\right]<N$. Taking advantage of the diagonal structure of $\mathit{\Sigma }\left(\alpha \right)$ and orthogonality of the matrices $\mathit{V}\left(\alpha \right)$ and $\mathit{U}\left(\alpha \right)$, it may be simply proved that the parameter ${\overline{\mathit{g}}}_K^{\lambda }\left(\alpha \right)$ (27) is given by

$${\overline{\mathit{g}}}_K^{\lambda }\left(\alpha \right)=\mathit{V}\left(\alpha \right){\mathit{\Lambda }}_{\lambda }\left(\alpha \right) \mathit{U}{\left(\alpha \right)}^T {\overline{\mathit{G}}}_N,$$

(48)

where $K\times N$ diagonal matrix ${\mathit{\Lambda }}_{\lambda }\left(\alpha \right)$ is as follows:

$${\mathit{\Lambda }}_{\lambda }\left(\alpha \right)=diag\left(\frac{{\sigma }_1\left(\alpha \right)}{{\left[{\sigma }_1\left(\alpha \right)\right]}^2+\lambda },\dots ,\frac{{\sigma }_{r\left(\alpha \right)}\left(\alpha \right)}{{\left[{\sigma }_{r\left(\alpha \right)}\left(\alpha \right)\right]}^2+\lambda },0,\dots ,0\right).$$

(49)

Using SVD (47) and introducing $N$ dimensional vector $\mathit{Y}\left(\alpha \right)=\mathit{U}{\left(\alpha \right)}^T {\overline{\mathit{G}}}_N$, the GCV function (29) is also expressed by a convenient analytical formula:

$$V_{GCV}\left(\lambda ,\alpha \right)=\frac{\left[{\sum }_{i=1}^{r\left(\alpha \right)}\frac{{\lambda }^2{\left[y_i\left(\alpha \right)\right]}^2}{{\left({\left[{\sigma }_i\left(\alpha \right)\right]}^2+\lambda \right)}^2}+{\sum }_{i=r\left(\alpha \right)+1}^N {\left[y_i\left(\alpha \right)\right]}^2\right]}{{\left[N-r\left(\alpha \right)+{\sum }_{i=1}^{r\left(\alpha \right)}\frac{\lambda }{\left({\left[{\sigma }_i\left(\alpha \right)\right]}^2+\lambda \right)}\right]}^2},$$

(50)

as a function of the singular values ${\sigma }_i\left(\alpha \right)$ and elements $y_i\left(\alpha \right)$ of the vector $\mathit{Y}\left(\alpha \right)$.

Similarly, the identification index $Q_{Nopt}\left(\alpha \right)$ (36) minimized in the upper level of the scheme is expressed using the SVD (47). Since, by virtue of (47) and (32) we have

$${\left(\mathit{\Psi }\left(\alpha \right)+{\lambda }_{GCV}\left(\alpha \right){\mathbb{I}}_{N,N}\right)}^{-2}=\mathit{U}\left(\alpha \right){\left[\mathit{\Sigma }\left(\alpha \right) \mathit{\Sigma }{\left(\alpha \right)}^T+{\lambda }_{GCV}\left(\alpha \right){\mathbb{I}}_{N,N}\right]}^{-2}\mathit{U}{\left(\alpha \right)}^T,$$

the index ${\widehat{Q}}_N\left(\alpha \right)$ is expressed as

$$Q_{Nopt}\left(\alpha \right)={\lambda }_{GCV}^2\left(\alpha \right)\left[{\sum }_{i=1}^{r\left(\alpha \right)}\frac{{\left[y_i\left(\alpha \right)\right]}^2}{{\left({\left[{\sigma }_i\left(\alpha \right)\right]}^2+{\lambda }_{GCV}\left(\alpha \right)\right)}^2}+{\sum }_{i=r\left(\alpha \right)+1}^N {\left[y_i\left(\alpha \right)\right]}^2\right].$$

(51)

3.5. Computational Algorithm for Model Identification

To design a numerical algorithm of the scheme, the communication between the levels should also be resolved. The computations must be arranged hierarchically in a two-level structure, i.e., for each iteration of the minimization procedure at the upper level, the whole numerical procedure must be realized for the lower-level GCV task (28). The complete computational procedure for determining the optimal model is given below.

• Step 0: Perform the experiment—stress relaxation test [1,2,3,28]—and record the measurements $\overline{G}\left(t_i\right)$, $i=1,\dots ,N$, of the relaxation modulus at times $t_i\ge 0$.
• Step 1: Determine the optimal regularization parameter ${\lambda }_{GCV}\left({\alpha }_{opt}\right)$ in the following two-level computations.
  • Step 1.0: Choose the initial point ${\alpha }^0$ for the numerical procedure applied to solve the upper-level task (37).
  • Step 1.1: Let ${\alpha }^m$ be the $m$-th iterate in the numerical procedure chosen to solve the upper-level task (37). For $\alpha ={\alpha }^m$, solve the lower-level minimization task (28) according to the chosen numerical optimization procedure and determine the regularization parameter ${\lambda }_{GCV}\left({\alpha }^m\right)$. The algebraic formula $V_{GCV}\left(\lambda ,\alpha \right)$ (50) is applied.
  • Step 1.2: Using ${\lambda }_{GCV}\left({\alpha }^m\right)$, compute, according to the numerical procedure selected to solve the upper-level task (37), with the index $Q_{Nopt}\left(\alpha \right)$ described by (51), the new parameter ${\alpha }^{m+1}$, which is the next approximation of ${\alpha }_{opt}$. If for ${\alpha }^{m+1}$ the stopping rule of the chosen numerical procedure is satisfied, i.e.,

    $$\left|{\alpha }^{m+1}-{\alpha }^m\right|\le {\epsilon }_1$$

    or

    $$\left|Q_{Nopt}\left({\alpha }^{m+1}\right)-Q_{Nopt}\left({\alpha }^m\right)\right|\le {\epsilon }_2,$$

    where ${\epsilon }_1$ and ${\epsilon }_2$ are preselected small positives, put ${\alpha }_{opt}={\alpha }^m$ as the optimal time-scale factor, ${\lambda }_{GCV}\left({\alpha }^m\right)$ as ${\lambda }_{GCV}\left({\alpha }_{opt}\right)$, and go to Step 2. Otherwise, return to Step 1.1 and continue the computations for $\alpha ={\alpha }^{m+1}$.
• Step 2: Compute the vector of the optimal model parameters ${\widehat{\mathit{g}}}_K$ according to (45) and the best model of the relaxation spectrum $H_{Kopt}\left(\tau \right)$ given by (46).

Remark 1.

The appealing feature of the scheme is that only the values of ${\lambda }_{GCV}\left({\alpha }^m\right)$, not the related parameters vector ${\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }^m\right)$ (34), are used for ${\alpha }^m$ in successive iterations of the numerical procedure solving the upper-level task (37).

Remark 2.

The regularization parameter ${\lambda }_{GCV}\left({\alpha }^m\right)$ resulting from the lower-level minimization task (28) in each iteration of the upper level is the solution of the GCV problem (28) for current ${\alpha }^m$. Thus, the respective vector ${\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }^m\right)$ (34) can be treated as an approximate solution of the overall identification problem.

Remark 3.

The selection of the initial time-scaling factor ${\alpha }^0$ in Step 1.0 may be based on the data concerning model applicability summarized in Table 1 and Table A1 or ${\alpha }^0$ can be selected by comparison, for different values of $\alpha$, the first basis function ${\varphi }_k\left(t,\alpha \right)$ (9), with the experiment results $\left\{\overline{G}\left(t_i\right)\right\}$.

Remark 4.

The SVD (47) of the matrix ${\mathit{\Phi }}_{N,K}\left(\alpha \right)$, of computational complexity $\mathcal{O}\left(N{K}^2\right)$ [57], must be computed only once for successive ${\alpha }^m$ generated in Step 1.1. SVD is accessible in the form of optimized numerical procedures in most commonly used computational packets.

The above procedure and communication between the levels are illustrated in Figure 3.

3.6. Analysis

Recently, a class of robust algorithms of the approximation of the continuous spectrum of relaxation frequencies by finite series of orthonormal functions has been developed and analyzed in detail in [40]. However, these results were derived using the orthogonality of the basis functions, so they are not applicable here. Therefore, both the analysis of the smoothing of the model and model accuracy for noisy measurements of the relaxation modulus must be carried out anew here.

3.6.1. Smoothness

The purpose of regularization relies on the stabilization of the resulting vector ${\overline{\mathit{g}}}_K^{\lambda }\left(\alpha \right)$ (27). Since the basis functions $h_k\left(\tau ,\alpha \right)$ (3) and (4) are such that $h_k\left(\tau ,\alpha \right)\le 1$ for any arguments, the following inequality

$$\underset{\tau \ge 0}{max}\left|H_K\left(\tau ,\alpha \right)\right|\le {\sum }_{k=0}^{K-1}\left|g_k\right|$$

holds for an arbitrary time-scale factor, which means that the smoothing of the vector of model parameters results in the limitation of the respective relaxation spectrum. The mechanism of the vector ${\overline{\mathit{g}}}_K^{\lambda }\left(\alpha \right)$ (27) stabilization via Tikhonov regularization is explained in many papers, for example, [20,40,55]. The following rule holds: the greater the regularization parameter $\lambda$ is, the more highly bounded the fluctuations of the vector ${\overline{\mathit{g}}}_K^{\lambda }\left(\alpha \right)$ are; see [40].

The norm ${\Vert H_K\left(\tau ,\alpha \right)\Vert }_2$ is a measure of smoothing of the relaxation spectrum model, where ${\Vert ·\Vert }_2$ also means the square norm in $L^2(0,\infty )$. In Appendix A.7, the following proposition is proved:

Proposition 1.

For an arbitrary time-scale factor $\alpha$ and arbitrary vector of model parameters ${\mathit{g}}_K$, we have

$${\Vert H_K\left(\tau ,\alpha \right)\Vert }_2^2={\mathit{g}}_K^T\mathit{\Gamma }\left(\alpha \right){\mathit{g}}_K=\frac{1}{2\alpha }{\mathit{g}}_K^T{\mathit{\Gamma }}_1{\mathit{g}}_K,$$

(52)

where $K\times K$ symmetric real matrices, $\mathit{\Gamma }\left(\alpha \right)$ described by (A39) and ${\mathit{\Gamma }}_1$ given by

$${\mathit{\Gamma }}_1=\left[\begin{array}{cccccc}1& \frac{e}{2}& \cdots & {\left(\frac{e}{2}\right)}^j\frac{j!}{j^j}& \cdots & {\left(\frac{e}{2}\right)}^{K-1}\frac{\left(K-1\right)!}{{\left(K-1\right)}^{K-1}}\\ \frac{e}{2}& 2{\left(\frac{e}{2}\right)}^2& \cdots & {\left(\frac{e}{2}\right)}^{1+j}\frac{\left(1+j\right)!}{j^j}& \cdots & {\left(\frac{e}{2}\right)}^K\frac{\left(K\right)!}{{\left(K-1\right)}^{K-1}}\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ {\left(\frac{e}{2}\right)}^k\frac{k!}{k^k}& {\left(\frac{e}{2}\right)}^{k+1}\frac{\left(k+1\right)!}{k^k}& \cdots & {\left(\frac{e}{2}\right)}^{k+j}\frac{\left(k+j\right)!}{k^k{j}^j}& \cdots & {\left(\frac{e}{2}\right)}^{k+K-1}\frac{\left(k+K-1\right)!}{k^k{\left(K-1\right)}^{K-1}}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ {\left(\frac{e}{2}\right)}^{K-1}\frac{(K-1)!}{{(K-1)}^{K-1}}& {\left(\frac{e}{2}\right)}^{K-1+j}\frac{\left(K\right)!}{{(K-1)}^{K-1}}& \cdots & {\left(\frac{e}{2}\right)}^{K-1+j}\frac{\left(K-1+j\right)!}{{(K-1)}^{K-1}{j}^j}& \cdots & {\left(\frac{e}{2}\right)}^{2K-2}\frac{\left(2K-2\right)!}{{(K-1)}^{2(K-1)}}\end{array}\right]$$

(53)

are a positive definite; matrix $\mathit{\Gamma }\left(\alpha \right)$ for an arbitrary time-scale factor $\alpha >0$.

The matrix ${\mathit{\Gamma }}_1$ is independent of the time-scale factor; only the multiplier $\frac{1}{2\alpha }$ in the last expression of (52) depends on $\alpha$. Since for any symmetric non-negative definite matrix, the eigenvalues and singular values are identical, by virtue of the Rayleigh–Ritz inequalities [58] (Lemma I):

$${\lambda }_{min}\left(\mathit{X}\right){\mathit{x}}^T\mathit{x}\le {\mathit{x}}^T\mathit{X}\mathit{x}\le {\lambda }_{max}\left(\mathit{X}\right){\mathit{x}}^T\mathit{x},$$

(54)

which hold for any $\mathit{x}ϵ{ℛ}^m$ and any symmetric matrix $\mathit{X}={\mathit{X}}^Tϵ{ℛ}^{m,m}$, where ${\lambda }_{min}\left(\mathit{X}\right)$ and ${\lambda }_{max}\left(\mathit{X}\right)$ are minimal and maximal eigenvalues of the matrix $\mathit{X}$, Equation (52) implies the following estimations:

$$\frac{1}{2\alpha }{\sigma }_{min}\left({\mathit{\Gamma }}_1\right){\Vert {\mathit{g}}_K\Vert }_2^2\le {\Vert H_K\left(\tau ,\alpha \right)\Vert }_2^2\le \frac{1}{2\alpha }{\sigma }_1\left({\mathit{\Gamma }}_1\right){\Vert {\mathit{g}}_K\Vert }_2^2,$$

where ${\sigma }_1\left({\mathit{\Gamma }}_1\right)$ and ${\sigma }_{min}\left({\mathit{\Gamma }}_1\right)$ denote the largest and the minimal singular values of matrix ${\mathit{\Gamma }}_1$ (53). Thus, the next result is derived.

Proposition 2.

For an arbitrary time-scale factor $\alpha$ and arbitrary vector of model parameters ${\mathit{g}}_K$, the following inequalities hold:

$$\frac{1}{\sqrt{2\alpha }}\sqrt{{\sigma }_{min}\left({\mathit{\Gamma }}_1\right)}{\Vert {\mathit{g}}_K\Vert }_2\le {\Vert H_K\left(\tau ,\alpha \right)\Vert }_2\le \frac{1}{\sqrt{2\alpha }}\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}{\Vert {\mathit{g}}_K\Vert }_2.$$

(55)

The square roots of the singular values ${\sigma }_1\left({\mathit{\Gamma }}_1\right)$ and ${\sigma }_{min}\left({\mathit{\Gamma }}_1\right)$ for $K=5,6,\dots 12$ are summarized in

Table 2. The square roots of the largest ${\sigma }_1\left({\mathit{\Gamma }}_1\right)$ and minimal ${\sigma }_{min}\left({\mathit{\Gamma }}_1\right)$ singular values of the matrix ${\mathit{\Gamma }}_1$ (53) for $K=4,5,\dots 12$ model summands.

K	4	5	6	7	8	9	10	11	12
$\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}$	3.666396	4.325185	4.927061	5.485792	6.010247	6.506534	6.979079	7.431223	7.865567
$\sqrt{{\sigma }_{min}\left({\mathit{\Gamma }}_1\right)}$	0.206481	0.087523	0.034642	0.013243	0.004953	0.001824	6.6358 × 10−4	2.3925 × 10−4	8.5615 × 10−5

. However, the lower bound of this norm is useful only for small $K$ and small time-scale factors.

Since $\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}$ grows with $K$, the greater the number of model summands there are, the greater the time-scaling factor should be to achieve pre-assumed multiplier $\frac{1}{\sqrt{2\alpha }}\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}$ in the estimation (55). In [59], a decreasing sequence of upper bounds for the largest singular value of a non-negative definite square matrix is constructed, given by Equation (19) in [59]. This result applied to the $K\times K$ matrix ${\mathit{\Gamma }}_1$ means that

$${\psi }_n\left({\mathit{\Gamma }}_1\right)=\frac{tr\left({\mathit{\Gamma }}_1\right)}{K}+{\left[\frac{{\left(K-1\right)}^{2n-1}}{{\left(K-1\right)}^{2n-1}+1}tr\left[{\left({\mathit{\Gamma }}_1-\frac{tr\left({\mathit{\Gamma }}_1\right)}{K}{\mathbb{I}}_{K,K}\right)}^{2n}\right]\right]}^{\frac{1}{2n}},$$

(56)

is a decreasing sequence of upper bounds for ${\sigma }_1\left({\mathit{\Gamma }}_1\right)$. The right inequality in (55) can be weakened to the following:

Proposition 3.

For an arbitrary time-scale factor $\alpha$ and arbitrary vector of model parameters ${\mathit{g}}_K$, the sequence of inequalities hold:

$${\Vert H_K\left(\tau ,\alpha \right)\Vert }_2\le \frac{1}{\sqrt{2\alpha }}\sqrt{{\psi }_n\left({\mathit{\Gamma }}_1\right)}{\Vert {\mathit{g}}_K\Vert }_2$$

(57)

for $n=1,2,\dots$, with the coefficients ${\psi }_n\left({\mathit{\Gamma }}_1\right)$ (56).

In

Table A2. The square roots of the upper bounds ${\psi }_n\left({\mathit{\Gamma }}_1\right)$ (56) and the square roots of the largest ${\sigma }_1\left({\mathit{\Gamma }}_1\right)$ and minimal ${\sigma }_{min}\left({\mathit{\Gamma }}_1\right)$ singular values of the matrix ${\mathit{\Gamma }}_1$ (53) for $K=4,5,\dots 12$ model summands.

$\mathit{K}$	4	5	6	7	8	9	10	11	12
$\mathit{n}$	$\sqrt{{\mathit{\psi }}_{\mathit{n}}\mathbf{\left(}{\mathit{\Gamma }}_{\mathbf{1}}\mathbf{\right)}}$
1	3.680199	4.353830	4.974923	5.556319	6.106199	6.630173	7.132299	7.615642	8.082583
2	3.671711	4.330308	4.931579	5.489724	6.013827	6.510086	6.982965	7.435811	7.871217
3	3.667383	4.325690	4.927344	5.485966	6.010365	6.506625	6.979169	7.431342	7.865747
4	3.666555	4.325232	4.927078	5.485799	6.010251	6.506536	6.979081	7.431227	7.865575
5	3.666421	4.325190	4.927062	5.485792	6.010247	6.506534	6.979079	7.431224	7.865568
6	3.666400	4.325186	4.927061	5.485792	6.010247	6.506534	6.979079	7.431223	7.865567
7	3.666397	4.325186	4.927061	5.485792	6.010247	6.506534	6.979079	7.431223	7.865567
8	3.666396	4.325186	4.927061	5.485792	6.010247	6.506534	6.979079	7.431223	7.865567
9	3.666396	4.325186	4.927061	5.485792	6.010247	6.506534	6.979079	7.431223	7.865567
10	3.666396	4.325186	4.927061	5.485792	6.010247	6.506534	6.979079	7.431223	7.865567
$\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}$	3.666396	4.325185	4.927061	5.485792	6.010247	6.506534	6.979079	7.431223	7.865567

in the Appendix B, the sequences $\sqrt{{\psi }_n\left({\mathit{\Gamma }}_1\right)}$ are summarized for $n=1,2,\dots ,10$ and $K=4,5,\dots 12$ model summands; for comparison, in the last row, $\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}$ are given. The sequence $\sqrt{{\psi }_n\left({\mathit{\Gamma }}_1\right)}$ quickly decreases to $\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}$; already, the sixth to eighth estimates equal $\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}$. However, it should be remembered that the right inequality in (55) and (57) only give the upper bounds of the norm ${\Vert H_K\left(\tau ,\alpha \right)\Vert }_2$.

To summarize, the smoothness of the optimal solution ${\overline{\mathit{g}}}_K^{\lambda }\left(\alpha \right)$ (27) of discrete regularized problem (26) guarantees that the fluctuations in the respective spectrum of relaxation, in particular the resulting spectrum of relaxation ${\widehat{H}}_K\left(\tau \right)$ (46), are also bounded. The time-scale factor $\alpha$ also affects the smoothness of the spectrum model.

3.6.2. Convergence and Noise Robustness

The relaxation spectrum recovery from experimental data is an inverse ill-posed problem, in which the identification index refers to the measured relaxation modulus, but not directly to the unknown relaxation spectrum $H\left(\tau \right)$. Therefore, we cannot estimate the model error ${\Vert H\left(\tau \right)-H_K\left(\tau ,\alpha \right)\Vert }_2$ directly. As a reference point for the model $H_K\left(\tau ,\alpha \right)$ (6) and (27), we will consider the model of the spectrum that we would obtain for the same time-scale factor $\alpha$ and regularization parameter $\lambda$ on the basis of ideal (undisturbed) measurements of the relaxation modulus:

$${\tilde{H}}_K\left(\tau ,\alpha \right)={\sum }_{k=0}^{K-1}{g}_K^{\lambda }\left(\alpha \right)h_k\left(\tau ,\alpha \right),$$

(58)

where ${\mathit{g}}_K^{\lambda }\left(\alpha \right)$ is the vector of the regularized solution of (26) given by (compare (27))

$${\mathit{g}}_K^{\lambda }\left(\alpha \right)={\left({\mathit{\Phi }}_{N,K}^T\left(\alpha \right){\mathit{\Phi }}_{N,K}\left(\alpha \right)+\lambda {\mathbb{I}}_{K,K}\right)}^{-1}{\mathit{\Phi }}_{N,K}^T\left(\alpha \right){\mathit{G}}_N$$

(59)

for noise-free measurements of relaxation modulus ${\mathit{G}}_N={\left[\begin{array}{ccc}G\left(t_1\right)& \cdots & G\left(t_N\right)\end{array}\right]}^T$; c.f., Equation (24). In Appendix A.8, the following estimations are derived:

Proposition 4.

For an arbitrary time-scale factor $\alpha$ and arbitrary regularization parameter $\lambda$, the error between the relaxation time spectrum models $H_K\left(\tau ,\alpha \right)$ (6) and (27) and ${\tilde{H}}_K\left(\tau ,\alpha \right)$ (58) and (59) is estimated using the following inequality:

$${\Vert H_K\left(\tau ,\alpha \right)-{\tilde{H}}_K\left(\tau ,\alpha \right)\Vert }_2\le \frac{1}{\sqrt{2\alpha }}\sqrt{{\sigma }_1\left({\mathit{\Gamma }}_1\right)}{\left[{\sum }_{i=1}^{r\left(\alpha \right)}\frac{{\left[{\sigma }_i\left(\alpha \right)\right]}^4}{{\left\{{\left[{\sigma }_i\left(\alpha \right)\right]}^2+\lambda \right\}}^4}\right]}^{\frac{1}{4}}{\Vert {\mathit{z}}_N\Vert }_2$$

(60)

where ${\mathit{z}}_N={\left[\begin{array}{ccc}z\left(t_1\right)& \cdots & z\left(t_N\right)\end{array}\right]}^T$ is the vector of measurement noises.

According to inequality (60), the accuracy of the spectrum approximation depends both on the measurement noises and regularization parameter and on the singular values ${\sigma }_1\left(\alpha \right),\dots ,{\sigma }_{r\left(\alpha \right)}\left(\alpha \right)$ of the matrix ${\mathit{\Phi }}_{N,K}\left(\alpha \right)$ (47) depending on the time-scale factor. By (60), the spectrum $H_K\left(\tau ,\alpha \right)$ tends to the noise-free spectrum ${\tilde{H}}_K\left(\tau ,\alpha \right)$ in each relaxation time $\tau$, at which they are both continuous, linearly with respect to the norm ${\Vert {\mathit{z}}_N\Vert }_2$, as ${\Vert {\mathit{z}}_N\Vert }_2\to 0$.

3.7. Example

Consider a viscoelastic material of the relaxation spectrum described by the double-mode Gauss-like distribution:

$$H\left(\tau \right)=\left[{\beta }_1{e}^{-{\left(\frac{1}{\tau }-m_1\right)}^2/q_1}+{\beta }_2{e}^{-{\left(\frac{1}{\tau }-m_2\right)}^2/q_2}\right]/\tau ,$$

(61)

where the parameters are as follows: ${\beta }_1=467\mathrm{Pa}·\mathrm{s}$, $m_1=0.0037{\mathrm{s}}^{-1}$, $q_1=1.124261\times {10}^{-6}{\mathrm{s}}^{-2}$, ${\beta }_2=39\mathrm{Pa}·\mathrm{s}$, $m_2=0.045{\mathrm{s}}^{-1}$ and $q_2=1.173\times {10}^{-3}{\mathrm{s}}^{-2}$. Spectra of this type are tested at the stage of developing new identification methods, for example, in [23] (Figure 2), [24] (Figures 9, 11 and 17) and [25] (Figures 2, 3, 6, 7–11 and 14), because they describe the viscoelastic properties of various materials, in particular: polymers [60] (Figures 4b and 8b), polyacrylamide gels [28] (Figure A4), glass [61] (Figure 2), wood [13] and sugar beet [18] (Figure 2). It is shown in Appendix A.9 that the corresponding ‘real’ relaxation modulus is

$$\begin{gathered} G\left(t\right)=\frac{\sqrt{\pi }}{2}[{\beta }_1\sqrt{q_1} e^{\frac{1}{4}{t}^2{q}_1-m_{1}t}erfc\left(\frac{\frac{1}{2}t{q}_1-m_1}{\sqrt{q_1}}\right)+ \\ {\beta }_2\sqrt{q_2} e^{\frac{1}{4}{t}^2{q}_2-m_{2}t}erfc\left(\frac{\frac{1}{2}t{q}_2-m_2}{\sqrt{q_2}}\right)], \end{gathered}$$

(62)

where the complementary error function $erfc\left(x\right)$ is defined by [62]:

$$erfc\left(x\right)=\frac{2}{\sqrt{\pi } }{\int }_x^{\infty }{e}^{-z^2}dz.$$

(63)

In the experiment, $N=5000$ sampling instants $t_i$ were generated with the constant period in the time interval $\mathcal{T}=\left[0,1550\right]$ seconds selected in view of the course of the modulus $G\left(t\right)$ (62). Additive measurement noises $z\left(t_i\right)$ were selected independently by random choice with uniform distribution on the interval $\left[-0.005, 0.005\right] \mathrm{Pa}$. The ‘real’ spectrum (61), modulus (62) and the basis functions $h_k\left(\tau ,\alpha \right)$, ${\varphi }_k\left(t,\alpha \right)$ were simulated in Matlab R2022a using the special functions besselk and erfc. For the singular value decomposition procedure, svd was applied. New calculation algorithms of the modified Bessel function of the second kind are constantly being developed; recently, an algorithm for parallel calculation was proposed [63].

3.7.1. Optimal Models

For $K=3,4,\dots ,12$, the optimal time-scaling factors ${\alpha }_{opt}$ were determined via the proposed two-level identification scheme and are given in

Table 3. The parameters of the optimal models in the example: optimal time-scale factors ${\alpha }_{opt}$, upper ${\alpha }_{min}$ and lower ${\alpha }_{max}$ bounds of the interval $\left[{\alpha }_{min},{\alpha }_{max}\right]$ of suboptimal time-scale factors defined by inequality (66) for ${\epsilon }_{\alpha }=0.1$, regularization parameters ${\lambda }_{GCV}\left({\alpha }_{opt}\right)$, optimal identification indices $Q_{Nopt}\left({\alpha }_{opt}\right)$ defined in (37), square norms: ${\Vert {\widehat{\mathit{g}}}_K\Vert }_2$ of the vector of optimal model parameter ${\widehat{\mathit{g}}}_K$ (45) and ${\Vert H_{Kopt}\left(\tau \right)\Vert }_2$ of the optimal relaxation spectrum model $H_{Kopt}\left(\tau \right)$ (46), relaxation spectrum approximation error $ER1\left(K\right)$ defined in (64) and the distance between successive approximations $H_{Kopt}\left(\tau \right)$ measured by $ER2\left(K\right)$ (65).

$\mathit{K}$	${\mathit{\alpha }}_{\mathit{o}\mathit{p}\mathit{t}}$ $\left[{\mathbf{s}}^{-1}\right]$	${\mathit{\alpha }}_{\mathit{m}\mathit{i}\mathit{n}}$ $\left[{\mathbf{s}}^{-1}\right]$	${\mathit{\alpha }}_{\mathit{m}\mathit{a}\mathit{x}}$ $\left[{\mathbf{s}}^{-1}\right]$	${\mathit{\lambda }}_{\mathit{G}\mathit{C}\mathit{V}}\left({\mathit{\alpha }}_{\mathit{o}\mathit{p}\mathit{t}}\right)$ $[-]$	${\mathit{Q}}_{\mathit{N}\mathit{o}\mathit{p}\mathit{t}}\left({\mathit{\alpha }}_{\mathit{o}\mathit{p}\mathit{t}}\right)$ $\left[{\mathbf{Pa}}^2\right]$	${\mathbf{\Vert }{\widehat{\mathit{g}}}_{\mathit{K}}\mathbf{\Vert }}_2$ $\left[\mathbf{Pa}\right]$	${\mathbf{\Vert }{\mathit{H}}_{\mathit{K}\mathit{o}\mathit{p}\mathit{t}}\left(\mathit{\tau }\right)\mathbf{\Vert }}_2$ $\left[\mathbf{Pa}\right]$	$\mathit{E}\mathit{R}1\left(\mathit{K}\right)$ $\left[\mathbf{\%}\right]$	$\mathit{E}\mathit{R}2\left(\mathit{K}\right)$ $\left[{\mathrm{Pa}}^2·\mathbf{s}\right]$
3	0.00520	0.00205	0.0080	0.0445	8.63505 × 10−4	0.7055	13.0248	90.459	12.876
4	0.01675	0.0164	0.01745	0.0063	3.43945 × 10−5	7.3485	16.1607	39.160	3.1049
5	0.02025	0.0192	0.0213	0.0071	2.71552 × 10−5	4.6724	17.1465	35.638	2.4592
6	0.02375	0.0220	0.0255	0.0078	2.48511 × 10−5	3.6493	18.0443	33.986	0.7563
7	0.02655	0.0234	0.0290	0.0083	2.48256 × 10−5	2.8846	18.5088	33.364	0.1432
8	0.02865	0.0234	0.0318	0.0089	2.51617 × 10−5	2.3555	18.5679	32.824	0.1692
9	0.03005	0.0241	0.03425	0.0099	2.52412 × 10−5	2.0639	18.3756	32.701	0.0111
10	0.03215	0.0255	0.0367	0.0109	2.51143 × 10−5	1.9058	18.4224	32.631	0.0336
11	0.03390	0.0276	0.03915	0.0122	2.48521 × 10−5	1.8020	18.3498	32.879	0.0327
12	0.03670	0.02935	0.04195	0.0127	2.44452 × 10−5	1.7198	18.4432	32.919	0.0328

together with the related regularization parameters ${\lambda }_{GCV}\left({\alpha }_{opt}\right)$. Next, the vectors of optimal model parameters ${\widehat{\mathit{g}}}_K={\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }_{opt}\right)$ (45) were computed and are given in

Table A3. Optimal parameters ${\widehat{\mathit{g}}}_K={\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left(\widehat{\alpha }\right)$ (45) of the relaxation spectrum models from Example 1; the elements of the vectors ${\widehat{\mathit{g}}}_K$ are expressed in $\left[\mathrm{Pa}\right]$.

${\widehat{\mathit{g}}}_{\mathit{K}}$
$\mathit{K}=3$	$\mathit{K}=4$	$\mathit{K}=5$	$\mathit{K}=6$	$\mathit{K}=7$	$\mathit{K}=8$	$\mathit{K}=9$	$\mathit{K}=10$	$\mathit{K}=11$	$\mathit{K}=12$
0.50819	−0.32787	−0.31484	−0.28238	−0.26421	−0.26104	−0.26344	−0.25958	−0.25722	−0.24893
0.48916	3.25406	2.44705	1.86341	1.57535	1.44701	1.38846	1.26949	1.18396	1.042852
−0.01177	−5.37283	−1.74791	−0.21483	0.21735	0.33298	0.38151	0.53024	0.62627	0.77959
	3.79955	−1.90531	−1.81625	−1.22415	−0.94252	−0.85336	−0.72742	−0.63340	−0.49513
		3.00974	−0.42686	−0.80274	−0.65867	−0.57137	−0.58067	−0.58715	−0.58709
			2.49772	0.35990	0.00038	0.018141	−0.06714	−0.11806	−0.19722
				1.85736	0.61157	0.39343	0.27044	0.22125	0.12879
					1.25669	0.58048	0.40004	0.34372	0.26284
						0.75498	0.46114	0.34725	0.26954
							0.59592	0.35547	0.25643
								0.46177	0.31367
									0.49927

in Appendix B; the elements of these vectors are both negative and positive. The square norms ${\Vert {\widehat{\mathit{g}}}_K\Vert }_2$ and ${\Vert H_{Kopt}\left(\tau \right)\Vert }_2$ are also enclosed in Table 3 as the measures of the solution smoothness. The norm ${\Vert H_{Kopt}\left(\tau \right)\Vert }_2$ of the optimal model $H_{Kopt}\left(\tau \right)$ (46) was determined through (52) based on ${\widehat{\mathit{g}}}_K$. For the ‘real’ spectrum $H\left(\tau \right)$ (61), the norm ${\Vert H\left(\tau \right)\Vert }_2=19.2562\mathrm{Pa}$. The approximation error between $H\left(\tau \right)$ (61) and their models $H_{Kopt}\left(\tau \right)$ (46) was estimated via relative integral error $ER1\left(K\right)$, defined by:

$$ER1{\left(K\right)}^2=\frac{{\Vert H\left(\tau \right)-H_{Kopt}\left(\tau \right)\Vert }_2^2}{{\Vert H\left(\tau \right)\Vert }_2^2}=\frac{{\int }_0^{\infty }{\left[H\left(\tau \right)-H_{Kopt}\left(\tau \right)\right]}^{2}d\tau }{{\int }_0^{\infty }H{\left(\tau \right)}^{2}d\tau },$$

(64)

which is expressed in a percentage in the penultimate column of Table 3. For the bimodal spectrum, the values of this error of the order of 33% are not surprising in the context of the ill-conditioned inverse problem. To compare the approximations $H_{Kopt}\left(\tau \right)$ (46) obtained for successive integer $K$, the square index $ER2\left(K\right)$ defined by the distance:

$$ER2\left(K\right)={\Vert H_{\left(K+1\right)opt}\left(\tau \right)-H_{Kopt}\left(\tau \right)\Vert }_2^2={\int }_0^{\infty }{\left[H_{\left(K+1\right)opt}\left(\tau \right)-H_{Kopt}\left(\tau \right)\right]}^{2}d\tau ,$$

(65)

is applied; see the last column of Table 3. The optimal indices $Q_{Nopt}\left({\alpha }_{opt}\right)$ (37) are given in Table 3, too. The optimal models $H_{Kopt}\left(\tau \right)$ (46) and the ‘real’ spectrum $H\left(\tau \right)$ (61) are plotted in Figure 4 for $K=3,4,\dots 12$. The analysis of the data from Table 3 and, in particular, the analysis of Figure 4b–d show that increasing the number of model components above $K=8$ does not significantly affect the course of the model $H_{Kopt}\left(\tau \right)$ and its accuracy. This is emphasized by the values of the $ER1\left(K\right)$ (64) and $ER2\left(K\right)$ (65) indices; in particular, the integral square error $ER2\left(K\right)$ between successive $H_{Kopt}\left(\tau \right)$ and $H_{\left(K+1\right)opt}\left(\tau \right)$ spectrum approximations, which relates to the square norm ${\Vert H\left(\tau \right)\Vert }_2^2$ for $K\ge 8$ does not exceed 0.046%, and for $K\ge 10$, it is of the order of 0.009%.

In Figure 5, the optimal models of the relaxation modulus $G_K\left(t,{\alpha }_{opt}\right)$ computed for ${\widehat{\mathit{g}}}_K$ (45) according to (7) are plotted, where the measurements $\overline{G}\left(t_i\right)$ of the ‘real’ modulus $G\left(t\right)$ (62) are also marked. The optimal models $G_K\left(t,{\alpha }_{opt}\right)$ have been well fitted to the experimental data, as indicated by the mean-square model errors $Q_{Nopt}\left({\alpha }_{opt}\right)/N$, which for $8\le K\le 12$ vary in the range from $4.889·{10}^{-9}{\mathrm{Pa}}^2$ to $5.0482·{10}^{-9}{\mathrm{Pa}}^2$. Thus, models $G_K\left(t,{\alpha }_{opt}\right)$ for different $K$ practically coincide with the measurement points and with each other; see Figure 5.

3.7.2. Optimal and Sub-Optimal Time-Scale Factors

The identification index $Q_{Nopt}\left(\alpha \right)$ (35) minimized in the upper-level task (37) as a function of time-scale factor $\alpha$ is plotted in Figure 6a–d for $K=6,8,10,12$. The optimal parameters ${\alpha }_{opt}$ are marked. The analysis of these plots suggests that there is such a neighborhood of ${\alpha }_{opt}$, namely a closed interval $\left[{\alpha }_{min},{\alpha }_{max}\right]$, that ${\alpha }_{min}\le {\alpha }_{opt}\le {\alpha }_{max}$ and for each $\alpha \in \left[{\alpha }_{min},{\alpha }_{max}\right]$, the identification index $Q_{Nopt}\left(\alpha \right)$ differs from the minimal $Q_{Nopt}\left({\alpha }_{opt}\right)$ not more than ${\epsilon }_{\alpha }·Q_{Nopt}\left({\alpha }_{opt}\right)$, i.e.,

$$Q_{Nopt}\left(\alpha \right)-Q_{Nopt}\left({\alpha }_{opt}\right)\le {\epsilon }_{\alpha }·Q_{Nopt}\left({\alpha }_{opt}\right),$$

(66)

where ${\epsilon }_{\alpha }$ is a small positive number. Inequality (66) means the deterioration of the model error is not greater than ${\epsilon }_{\alpha }$ percent of the optimal $Q_{Nopt}\left({\alpha }_{opt}\right).$ Any parameter $\alpha$ from the interval $\left[{\alpha }_{min},{\alpha }_{max}\right]$ is a suboptimal time-scale factor. In Figure 6, ${\epsilon }_{\alpha }=0.1$, which means a 10% level of sub-optimality, is assumed, and the values of ${\alpha }_{min}$ and ${\alpha }_{max}$ are marked on the small subfigures. They are also given in Table 3. With the increase in the number of model components, the optimal parameter ${\alpha }_{opt}$ increases and the range of sub-optimal time-scale factors shifts. For an exemplary number of model components, $K=9$, the optimal model $H_{Kopt}\left(\tau \right)$ (46) and models $H_K\left(\tau ,\alpha \right)$ (6), optimal in the sense of lower-level task (28) for suboptimal parameters ${\alpha }_{min}$ and ${\alpha }_{max}$, are plotted in Figure 7a; the ‘real’ spectrum $H\left(\tau \right)$ (61) is presented, too. The corresponding optimal regularization parameters are as follows: ${\lambda }_{GCV}\left({\alpha }_{min}\right)=0.0167$ and ${\lambda }_{GCV}\left({\alpha }_{max}\right)=0.0094$. For respective vectors of optimal model parameters, we have ${\Vert {\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }_{min}\right)\Vert }_2=2.5114\mathrm{Pa}$ and ${\Vert {\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }_{max}\right)\Vert }_2=2.5119\mathrm{Pa}$. The norms of the relaxation spectra are ${\Vert H_K\left(\tau ,{\alpha }_{min}\right)\Vert }_2=17.4472\mathrm{Pa}$ and ${\Vert H_K\left(\tau ,{\alpha }_{max}\right)\Vert }_2=20.2069\mathrm{Pa}$. In Figure 7b, the ‘real’ modulus $G\left(t\right)$ (62) and the models: $G_K\left(t,{\alpha }_{opt}\right)$ (7) and (45) and $G_K\left(t,{\alpha }_{min}\right)$, $G_K\left(t,{\alpha }_{max}\right)$ computed for ${\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }_{min}\right)$, ${\widehat{\mathit{g}}}_K^{{\lambda }_{GCV}}\left({\alpha }_{max}\right)$ according to (7) are plotted. However, the relaxation moduli are almost identical (Figure 7b), the spectrum models differ (Figure 7a), which emphasizes the importance of the time-scale factor optimal selection. Too strong smoothing of the relaxation spectrum models, with a simultaneous very good fit to the experimental data of the relaxation modulus models, indicates the need to modify the quality index in the regularized task (26). For example, the square term ${\Vert {\mathit{g}}_K\Vert }_2^2$ in the second component of the objective function can be replaced by the quadratic form ${\mathit{g}}_K^T\mathit{W}{\mathit{g}}_K$, and a positive definite weight matrix $\mathit{W}$ or regularized weighted least-squares approach [41] can be applied. This will be the subject of further work.

4. Conclusions

In this paper, a new robust hierarchical algorithm for the identification of the relaxation time spectrum, which combines the technique of an expansion of a function into a series of basis functions with the least-squares regularized identification and the optimal choice of the time-scale factor, has been derived. The task of determining the best ‘regularized’ model was solved at the lower level, while the optimal time-scale factor was selected on the upper level of the identification scheme. The continuous spectrum of relaxation times was approximated by finite series of power–exponential basis functions, while the components of the relaxation modulus model were proven to be described by the product of power of time and the modified Bessel function of the second kind. In the present paper, the problem of the optimal choice of the time-scale factor to ensure the best fit of the model to experimental data has been formulated and solved for the first time in the context of the relaxation spectrum identification.

The necessary optimality conditions both for the optimal regularized least-squares identification task and the problem of the optimal selection of the time-scale factor were derived in the form of nonlinear algebraic equations. The main properties of the basis functions of the relaxation spectrum and modulus models, their positive definiteness, convenient upper bounds, monotonicity, asymptotic properties and wide range of applicability for different time-scale factors indicated the possibility of using the proposed model and identification algorithm to determine the spectrum of a wide class of viscoelastic materials.

The overly strong smoothing of the relaxation spectrum models in the example, with a very good fit to the experimental data of the relaxation modulus models, indicates the need to modify the quality index in the lower level identification task. An introduction of a weight matrix in the second component of the objective function or a direct application of regularized weighted least-squares should be investigated. Another solution is the selection of such a spectrum model which guarantees the assumed level of smoothing and optimal adjustment to the relaxation modulus measurement data. These approaches will be the subject of further work.

The presented scheme of the relaxation spectrum identification can be easily modified for retardation spectrum recovery from creep compliance measurements obtained in the standard creep test.