A Semi-Explicit Algorithm for Parameters Estimation in a Time-Fractional Dual-Phase-Lag Heat Conduction Model

: This paper presents a new semi-explicit algorithm for parameters estimation in a time-fractional generalization of a dual-phase-lag heat conduction model with Caputo fractional derivatives. It is shown that this model can be derived from a general linear constitutive relation for the heat transfer by conduction when the heat conduction relaxation kernel contains the Mittag– Leffler function. The model can be used to describe heat conduction phenomena in a material with power-law memory. The proposed algorithm of parameters estimation is based on the time integral characteristics method. The explicit representations of the thermal diffusivity and the fractional analogues of the thermal relaxation time and the thermal retardation are obtained via a Laplace transform of the temperature field and utilized in the algorithm. An implicit relation is derived for the order of fractional differentiation. In the algorithm, this relation is resolved numerically. An example illustrates the proposed technique.


Introduction
Fourier's law is a fundamental phenomenological relation in heat transfer theory that describes heat conduction in solids and fluids.However, it is not applicable for modelling the heat transfer in ultrafast processes (e.g., during laser heating and cooling [1]), on micro/nanoscales (e.g., heating of carbon nanotubes [2]), and in some complex media (e.g., heat conduction in biological tissues [3] and materials with a non-homogeneous inner structure [4]).For this reason, various linear and nonlinear generalizations of Fourier's law have been proposed by many researchers over the past two centuries.Historically, the first one is the Maxwell-Cattaneo-Vernotte (MCV) heat conduction model [5,6].It contains an additional differential term to describe thermal relaxation processes.This model leads to the hyperbolic heat conduction equation and overcomes the problem of infinite heat propagation speed in Fourier's law.A more common non-Fourier heat conduction model with a time delay in the heat flux is the single-phase-lag (SPL) model [7,8].Note that the MCV model can be considered as a linearisation of the SPL model with respect to the relaxation time.A further extension of the SPL model was developed by Tzou [9], who introduced an additional phase lag associated with the temperature gradient.This model is known as the dual-phase-lag (DPL) heat conduction model and permits one to take into account the microstructural effects during ultrafast heat transfer in a macroscopic formulation.Recently, time-fractional heat conduction has become a new branch of the heat transfer theory.Time-fractional models allow us to take into account memory thermal effects in a material by incorporating time derivatives of fractional orders [10,11] (usually, in the Caputo sense [11]) into the thermal constitutive relation.Povstenko [12] and Jiang and Xu [13] proposed a time-fractional generalization of Fourier's law and studied the corresponding fractional heat equation in different coordinate systems.Several timefractional generalizations of the MCV and SPL models were considered in [14][15][16], and time-fractional extensions of the DPL model were presented in [17,18].More detailed reviews of the non-Fourier heat conduction models can be found in [19][20][21][22][23][24][25].
Thermal conductivity is a unique property of a material in Fourier's law, whereas non-Fourier models usually include several material thermal properties.For example, the frequently used dual-phase-lag (DPL) heat conduction model incorporates the thermal conductivity, the thermal relaxation time, and the thermal retardation.As a result, the possibility of using a non-Fourier model in real-world applications is based on an ability to obtain the necessary thermal properties of a material.From a mathematical point of view, the problem of parameters estimation can be considered as an inverse coefficient problem, which is usually an ill-posed problem.Thus, the development of theoretical techniques for the estimation of thermal properties in non-Fourier models is a challenging problem of mathematical modelling in heat transfer theory.
There are two main paradigms in the field of thermal properties estimation, namely the deterministic and stochastic approaches [26].In the deterministic approach [27,28], the result of parameters estimation is obtained as a single solution but is not supposed to be exact.Usually, a deterministic technique reduces an inverse problem to an optimization problem.In view of measurement errors, the residual principle and regularization methods have to be used to obtain a stable solution [29][30][31].However, deterministic methods would be inappropriate when the number of parameters becomes quite large.In this case, the stochastic approach [32] is more preferable.In this approach, an inverse problem is considered in terms of probability in order to allow the uncertainty in measurements and models.The Bayesian framework is an effective stochastic technique for solving inverse problems in heat transfer [33].
The problems of parameters estimation for the MCV, SPL, and DPL models have been studied by various researchers.Since these models are linear with respect to all thermal parameters, the corresponding optimization problems arising during the solving of inverse problems are similar to classical ones for Fourier's model.Hence, the methods of parameters estimation developed for the classical heat conduction equation are also applicable for these models.For example, Kishore and Kumar [34] used the Levenberg-Marquardt algorithm for the parameters estimation of the DPL model.Mochnacki and Paruch [35] solved the problem of relaxation time identification in the MCV model by using evolutionary optimization algorithms.Liu and Lin [36] and Str ąkowska et al. [37] applied Laplace transform together with different optimization methods to identify the parameters of the DPL model.França and Orlande [38] studied the inverse parameter estimation problem for the DPL heat conduction model by using the Markov chain Monte Carlo method within the Bayesian framework.
In a time-fractional model, the orders of fractional differentiation are additional parameters that also have to be estimated.Since a time-fractional constitutive relation is nonlinear with respect to such orders, the problem of parameters estimation also becomes nonlinear.Ghazizadeh et al. [39], Yu et al. [40], and Mozafarifard et al. [41] solved the inverse time-fractional SPL heat conduction problem using the nonlinear parameter estimation technique based on the Levenberg-Marquardt method.Goudarzi et al. [42] solved a similar problem by employing the conjugate gradient inverse method.Sobhani et al. [43] considered the variable-order time fractional DPL model and estimated the two lagging times by using an inverse analysis based on the Levenberg-Marquardt algorithm.Qiao et al. [44] proposed modified hybrid Nelder-Mead simplex search and particle swarm optimization for parameters estimation in the time-fractional DPL model.Zheng et al. [45] used the Bayesian method to construct an algorithm to estimate the four parameters of the timefractional DPL model.Note that all of the mentioned approaches are based on a numerical solution of the corresponding direct problem.
The method of time integral characteristics (TICs) is an efficient technique for the estimation of constant parameters in linear models.It was proposed by Shatalov [46] at the end of the last century for solving the inverse coefficient problems of the heat conduction theory.The method is based on integral transformation of the initial-boundary value problem for the considered linear model on the time variable and solving the corresponding inverse coefficient problem in the transformed space.The absence of necessity to perform the inverse integral transform is a main advantage of this method.Also, it does not use a numerical solution of the direct problem.Later, this method was extended to classical and time-fractional diffusion models [47][48][49].
In this study, we focus on solving the problem of parameters estimation in a timefractional generalization of the DPL heat conduction model by the TIC method.The time-fractional dual-phase-lag (TFDPL) model was proposed by Xu and Jiang [17] to interpret the experiment results for processed meat.The Caputo time-fractional derivatives of two different orders are used in this model.The authors obtained an analytical solution for the corresponding bioheat transfer equation and solved the inverse problem for the estimation of model parameters by applying the nonlinear least-square method.The same model was used in [50] for treating the thermoelastic response of skin subjected to sudden temperature shock.In [51], a fundamental solution for the TFDPL heat conduction problem was obtained.Also, a TFDPL model with a single order of fractional differentiation was considered by several scholars.In [52], such a model was used for describing the heat conduction in a multi-layered spherical medium with azimuthal symmetry.In [18], a similar model with temperature jump boundary conditions was utilized for the numerical simulation of heat transfer in transistors.Numerical schemes for solving several TFDPL heat conduction problems was proposed in [53,54].
However, it is necessary to note that in all of the papers mentioned above the considered TFDPL models have been obtained from the classical DPL model by the formal replacing of the integer order time derivatives by their fractional analogues.In this paper, we overcome this weakness by proposing the derivation of the TFDPL model from a general linear constitutive relation for heat transfer by conduction.
The paper is organized as follows.Section 2 contains a brief description of the time integral characteristic method.Section 3 is devoted to the derivation of the TFDPL heat conduction model and corresponding non-Fourier heat conduction equation.A proposed semi-explicit algorithm for TFDPL model parameters estimation is described in Section 4.An illustrative example of using this algorithm is presented in Section 5. Section 6 contains the discussion of the obtained results and recommendations.Finally, Section 7 presents the concluding remarks.

Preliminaries
This section gives a brief description of the TIC method.This method was proposed for solving inverse problems of parameters estimation in linear evolution equations.Note that it is applicable only for constant coefficient equations.Also, it is assumed that the Laplace transform of the temperature field exists.
Let us illustrate the basic idea of the TIC method by a simple problem of the thermal diffusivity estimation.We consider the heat equation: Here, x and t are spatial and temporal variables, respectively, T(x, t) is the temperature field, and a is the thermal diffusivity.This equation is accompanied with the initial condition the boundary conditions and the additional internal condition Here, T 0 (t) and T l (t) are known functions.Then, the inverse coefficient problem is stated as follows: given the initial boundary value problem (1), (2), (3) and the additional condition (4), find the constant thermal diffusivity, a.
The method of TIC is based on an integral transformation of the temperature field with respect to time.The Laplace transform can be efficiently used for this purpose.The function T * L (p) = T * (L, p) is referred to as a time integral characteristic of the temperature field T at the point x = L.
The solution of ( 6) and ( 7) is Then, by using (8), we find the following explicit representation of the thermal diffusivity via TICs of the temperature field: A main advantage of the described technique is that there is no necessity to perform the inverse Laplace transform.If the temperature functions T 0 (t) and T l (t) are known exactly, the representation (9) gives the exact value of a for all permitted values of p.However, in practice such functions are usually measured in an experiment with some errors.Then, the Laplace parameter p should be considered as a regularization parameter and its value should be chosen in agreement with experimental errors.The explicit representation (9) permits one to obtain a linear estimate of the relative error for the thermal diffusivity as a function of p (see [46,48,49] for more details).The minimum of this function gives the optimal value of p.

The Time-Fractional Dual-Phase-Lag Heat Equation
A TFDPL heat conduction model can be obtained similarly to the time-fractional Zener model in the theory of linear fractional viscoelasticity [55].
To derive the model, we make the following assumptions.

Assumption 1.
A heat transfer medium is homogeneous and isotropic.
Assumption 2. The temperature history affects the heat flux (the thermal memory exists).

Assumption 3.
There is a linear relation between the heat flux and the temperature gradient (thermal relaxation process is linear).
Assumption 4. The heat conduction relaxation function r(t) is a differentiable, decreasing function of time.
Assumption 5.The relaxation function r(t) has a power-law decay in time (the so-called "heavy tail").
Assumption 6.The temperature gradient ∇T is a continuously differentiable function with respect to time and space variables.
Under Assumptions 1-3, a general linear constitutive relation for heat transfer by conduction in a medium is defined mathematically using a Riemann-Stieltjes integral as Here, q(t) ≡ q(x, t) is the heat flux vector, g(t) ≡ g(x, t) = −∇T is the temperature gradient, and r(t) is the heat conduction relaxation function, which does not depend on the spatial coordinate x.
In accordance with the physical principle of causality, the relaxation function r(t) is zero for negative time.Hence, the constitutive relation (10) takes the form It is easy to see that this equation reduces to Fourier's law, q = −k∇T if r(t) = k, where k is the thermal conductivity, which is a constant in time.Here, we additionally used the natural condition g(−∞) = 0, which means that the medium has finite total thermal energy over all times.
Based on the above Assumptions 4 and 6, we can rewrite (11) in a more convenient form: Let us now consider the case when the heat conduction relaxation function has the form where r 0 , r 1 , τ q are constants and is the Mittag-Leffler function (see, e.g., [11]).Note that the function r(t) defined by ( 13) has a power-law decay and therefore Assumption 5 is fulfilled.In the limiting case of α = 1, we have E 1 (z) = e z and relation (13) gives the heat conduction relaxation function for the classical DPL model.The Mittag-Leffler function ( 14) has the known property i.e., it is invariant with respect to the left-sided Caputo fractional derivative of order α.This fractional derivative reads Here, a I 1−α t is the left-sided fractional integral operator of order 1 − α.Letting a = −∞ in (16), we obtain the Caputo fractional operator C −∞ D α t for the whole axis R. Applying this operator to both sides of (12), we get Note that the existence of integrals in this expression is a consequence of Assumptions 5 and 6.Changing the order of integration in the last term of the above expression, we obtain Substituting r(t) given by ( 13) into (17), and using (15), we have On the other hand, relation ( 12) with ( 13) takes the form It is easy to see that the integrals in the last terms of ( 18) and ( 19) coincide and therefore can be excluded.As a result, we obtain the time-fractional constitutive relation Using the definition of the function g, this relation can be written as Here, k = r 0 is the thermal conductivity, τ q is the fractional analogue of the thermal relaxation time (the so-called phase lag in the heat flux), and τ T = τ q (1 + r 1 /r 0 ) is the fractional analogue of the thermal retardation (the so-called temperature gradient phase lag).
The constitutive relation (20) describes the heat conduction in a medium with full power-law memory.This relation is invariant with respect to translation in time, and therefore the time origin can be arbitrarily chosen.
Let us now assume additionally that there is no heat transfer in a medium for time t < 0.Then, relation (20) reduces to Note that this relation is not invariant with respect to translation in time, and the time origin is fixed in this case.Now we can obtain the TFDPL heat equation.The energy conservation law for a constant property material without heat sources can be written as where c is the specific heat and ρ is the density of the material.Combining ( 21) and ( 22), after simple algebra, we get where a = k/(cρ) is the thermal diffusivity.TFDPL heat Equation ( 22) models heat conduction in a medium with power-law memory and constant temperature field for time t < 0. Note that under the given assumptions we obtain a linear one-temperature heat conduction model with constant parameters and a single order of fractional differentiation, α.
A final brief remark should be made about the physical interpretation of the term "power-law memory".In complex heterogeneous media, a complexity of the heat transfer process can be associated with memory effects at the macroscopic description level.If these effects demonstrate power-law behaviour at large times, it is said that the medium or material has a power-law memory [56].It is worth noting that memory effects in such a medium are usually the results of averaging microscopic heterogeneities.In other words, a complex heterogeneous medium microscopically obeys Fourier's law, whereas macroscopically it is modelled as a homogeneous medium that obeys the non-Fourier's law with a power-law memory.Time-fractional derivatives are natural to describe such types of memory.Time-fractional heat conduction models were proposed, for example, for the porous media [43], for a carbon-carbon composite medium [57], for metal-oxidesemiconductor field-effect transistors [38], and for processed meat [17].

An Algorithm of TFDPL Model Parameters Estimation
Let us consider a one-dimensional case of TFDPL heat Equation (23) in a half space, namely This is a time-fractional equation of order α + 1 ∈ (1, 2), and therefore two initial conditions are needed for its unique solvability.We take them in the form where T 0 is a constant initial temperature.We will also assume that ( 24) is accompanied by the boundary conditions and by the additional internal condition Here, T 0 (t), T l (t) are known functions.We will consider the following inverse problem: given the initial boundary value problem (24), ( 25), (26) and the additional condition (27), find the constants a, τ T , τ q , α.For solving this problem, the TIC method can be efficiently used.
For convenience, we introduce a new function, θ(x, t) = T(x, t) − T 0 .Then, the problem ( 24)-( 27) takes the form Here, the functions θ 0 (t) = T 0 (t) − T 0 and θ l (t) = T l (t) − T 0 are known.The initial-boundary problem ( 28), ( 29), (30) after Laplace transform can be written as and (31) gives Here, θ * (x, p) denotes the Laplace transform of θ(x, t), which is defined by In (32), prime denotes differentiation with respect to x.The solution of ( 32), ( 33) is Using the additional condition (34), we obtain the main equation for parameters estimation, which can be written as where Note that the function Φ is linear with respect to a, b, τ q , and nonlinear with respect to α.
As mentioned in the Preliminaries section, the problem of finding the Laplace parameter p arises if the functions T 0 (t) and T l (t) are not known exactly.In the TIC method, the Laplace parameter p is assumed to be real and positive.Therefore, it is natural to assume that this parameter belongs to a finite interval, [p min , p max ] (0 < p min < p max < ∞).A discussion of different approaches to the estimation of p min and p max can be found in [48,49].Then, the considered inverse problem can be reduced to a minimization problem: This is a classical problem of finding a minimum of four variables function, f .The physical constraints are a > 0, b > 0, τ q > 0, and α ∈ (0, 1).In general, this problem can be solved numerically by using different optimization software.However, explicit TIC representations for the desired parameters can be obtained in a special case when the order of fractional differentiation α is known.Let us consider (37) as the unconstrained minimization problem.It is obvious that function f defined by ( 37) is a quadratic function in three variables, a, b, and τ q .The necessary conditions for local optimality read These conditions provide the system of linear equations where , and Note that the matrix A is symmetric.Using Cramer's rule, the solution of ( 38) can be written in the explicit form where and A i is the matrix formed by replacing the i-th column of A by the column vector B. Thus, the explicit representations (39) permit one to obtain the values of a, b, and τ q for a given value of α.
The representations (39) can also be used in the case of unknown α.Then, we have where In (40), the parameters a, b, and τ q are the functions of α, which are defined by (39).We thus obtain a single equation for α.Equation ( 40) is nonlinear and quite complex.Therefore, it should be solved numerically.
As a result, we can state the following semi-explicit algorithm for parameters estimation in the TFDPL heat conduction model:

3.
The order of fractional differentiation α is obtained by numerically solving nonlinear Equation (40).

4.
The desired parameters a, b, and τ q are calculated by the explicit TIC representations given by (39).
It is necessary to note that the proposed algorithm is based on the unconstrained minimization problem.As a result, the order of fractional differentiation α, which is obtained as the solution of (40), does not necessarily belong to the interval [0, 1].Then, the constrained minimization problem mentioned above should be considered and solved numerically.Note that usually this situation arises when the initial functions θ 0 (t) and θ l (t) have quite large errors (usually more than 5 %).
The considered problem of parameters estimation belongs to the class of inverse coefficient problems.Hence, it is an ill-posed problem in most cases.In the proposed algorithm, the stabilization of the solution is achieved by integration with respect to the Laplace parameter p.However, numerical experiments show that the solution is stable only if p max /p min = O(10 k ) with k ≥ 1.If p max /p min = O(1), the determinant det(A) is close to zero, and the corresponding approximate solution is unstable.Also, det(A) → 0 with p max → 0. An additional regularization is needed in this case.For example, the the Tikhonov regularization method can be used for this purpose.By applying this method to the system (38), we get a regularized system where A T is the transpose matrix of matrix A and ε is the regularization parameter.This parameter be found from the discrepancy principle if an upper error bound of the experimental temperature measurements is known.However, in this case we dramatically increase the computational cost of the algorithm and therefore remove the main advantage of the TIC approach.Thus, the parameters p min and p max should be chosen so that det(A) is not close to zero.Calculations show that the value of det(A) highly depends on α.Hence, the values of p min and p max should also be chosen dependently on α.Note that for a given interval [p min , p max ], the numerical value of det(A) can be calculated for any α ∈ (0, 1).If α is initially unknown, an appropriate interval [p min , p max ] can be found iteratively.

An Example
To illustrate the above algorithm, let us consider Equation ( 28) with Hence, b = 8 and Three different values of fractional order are considered: α = 0.75, 0.5, and 0.25.Graphs of the function ψ(p) defined by (42) for these values of α are shown in Figure 1.Denote by T(x, t) and T(x, t), the exact and perturbed temperature fields, respectively.Let where ∆ T is the upper error bound.Then, it is easy to prove (see, e.g., [46]) that Thus, the error of T * (x, p) is increased as the parameter p is decreased.The same is valid for the function ψ(p).
To simulate experimental errors, we approximated the function ψ(p) on the interval [p min , p max ] by polynomials of different degrees with respect to a new dependent variable P = ln p.
We utilized the relative error as a metric of accuracy.For a quantity q, it is defined by where q and q are exact and perturbed values of q, respectively.All computations were performed in Maple 17.
To illustrate the fact that p min and p max should depends on α, we firstly consider the fixed interval p ∈ [0.01, 10].
The following approximations of ψ(p) were constructed: Graphs of relative errors for these functions are plotted in Figure 2. It can be seen that the maximum value of relative error is approximately equal to 6.5% for ψ 1 , 2% for ψ 2 , and 0.5% for ψ 3 .

Figure 2. Graphs of relative errors for functions (43).
Table 1 contains the results of parameters estimation by using the explicit expressions given in (39) for α = 0.75 and approximate functions ψ i from (43).Note that, in this case, det(A) = O(10 −1 ).It can be seen that the relative errors of a and b are of the same order of magnitude as the corresponding relative errors of functions ψ i .The relative errors of τ T and τ q are also in good agreement with the relative errors of ψ 3 and ψ 2 .However, in the case of quite large initial error when the function ψ 1 is used, the error level of τ T and τ q is highly increased.Table 2 contains the results of parameters estimation for unknown α.In this case, Equation ( 40) was solved for each ψ i (i = 1, 2, 3).As can be seen from the table, the accuracy of α identification is highly depends on the error of input data.The same is valid for the other parameters.Case 2: α = 0.5.
The following approximations of ψ(p) were obtained: Graphs of relative errors for these functions are plotted in Figure 3.The maximum value of relative error is approximately equal to 5% for ψ 1 , 2.2% for ψ 2 , and 0.5% for ψ 3 .The results of parameters estimation are given in Table 3 for α = 0.5, and in Table 4 for initially unknown α.In general, the results of this case are close to the previous one.However, the magnitudes of relative errors for estimated parameters are greater than those obtained previously.In this case, det(A) = O(10 −3 ).Case 3: α = 0.25.
In this case, the following approximations of ψ(p) were constructed: Graphs of relative errors for the functions in (45) are plotted in Figure 4.The maximum value of relative error is approximately equal to 10% for ψ 1 , 1% for ψ 2 , and 0.1% for ψ 3 .Table 5 contains the results of parameters estimation by using (39) for α = 0.25 and approximations (45).In this case, we found det(A) = O(10 −5 ) and therefore the estimation results are highly sensitive to errors of input data.It follows from the table that if we use the ψ 1 function, having a relative error of, at most, ±10%, the proposed algorithm does not permit us to identify τ q with the chosen values of p min and p max .However, we obtained reasonable values for all of the desired parameters for the function ψ 3 .This demonstrates the stability of the algorithm.To explain the obtained results, we first consider det(A) as a function of α.The graphs of this function for different intervals [p min , p max ] are given in Figures 5 and 6.It can be seen that decreasing α decreases the value of det(A).Therefore, the sensitivity of the estimated parameters to measurement errors increases with decreasing the order of fractional differentiation.The simulation results given in the above tables confirm this conclusion.
As it follows from Figure 5, the value of det(A) highly depends on the value of p max , so that bigger values of det(A) correspond to bigger values of p max .The value of det(A) also depends on p min .However, the comparison of Figures 5 and 6 shows that this dependence is weaker than the previous one.Therefore, the large value of det(A) cannot be used as a criterion for finding the value of p min , especially in the case of small values of α.An essential feature of the heat conduction processes described by the TFDPL model is that the corresponding temperature field has a power-law asymptotic behaviour at large times, known as a "heavy tail".This means that such a temperature field varies significantly over a larger time interval than a temperature field described by the classical Fourier model, and this interval increases with decreasing α.In view of the final value theorem of the Laplace transform lim one can conclude that a smaller value of p min should be used in the algorithm for a smaller value of α.Indeed, for the function ψ(p) defined by (42), we have ψ ∈ (0.025, 1) for p ∈ (0, ∞) and any α ∈ (0, 1).However, for a fixed interval, p ∈ [0.01, 10], we obtain ψ ∈ [0.037, 0.928] for α = 1, ψ ∈ [0.046, 0.803] for α = 0.75, ψ ∈ [0.062, 0.567]) for α = 0.5, and ψ ∈ [0.089, 0.301] for α = 0.25.Hence, for α = 1 (the case of the classical DPL model) more than 90% of available information about the heat transfer process will be used in the algorithm if p ∈ [0.01, 10].Therefore, this interval will be close to the optimal one in this case.For this reason, the interval [0.01, 10] has been taken at the beginning of this example.However, less than 25% of available information about the heat transfer process was used in the algorithm for α = 0.25.Moreover, we lost all information about the "heavy tail".This was the main reason for high errors in parameters estimation and failure of the algorithm to identify τ q for the ψ 1 function.
To confirm this conclusion, let us consider the case of α = 0.25 with p min = 10 −6 and p max = 100.Then, we have ψ ∈ [0.0288, 0.9997] and det(A) = O(10 −1 ).The following three approximations of ψ(p) were constructed: Graphs of relative errors for these functions are plotted in Figure 7.It can be seen that the maximum value of relative error is approximately equal to 10% for ψ 1 , 0.9% for ψ 2 , and 0.5% for ψ 3 .Table 7 presents the results of parameters estimation by using the explicit expressions given in (39) for α = 0.25 and approximate functions ψ i from (46).A comparison of the results given in Table 7 with those given in Table 5 shows that for the enlarged interval p ∈ [10 −6 , 100] the proposed algorithm provides much better results than for p ∈ [0.01, 10].Table 8 contains the results of parameters estimation for unknown α.It can be seen that for the function ψ 3 with an error level of approximately 0.5%, the proposed algorithm permits one to estimate all parameters with high accuracy.However, for the more perturbed functions ψ 2 and ψ 3 , the obtained parameters have quite large errors.

Discussion and Recommendations
Now, let us briefly discuss the advantages and limitations of the algorithm described in Section 4.
The algorithm is based on the exact analytical solution of the considered initialboundary value problem for the TFDPL heat conduction equation in the Laplace space.As a result, there is no necessity to solve this problem numerically.Also, it is not necessary to perform the inverse Laplace transform.Note that most modern algorithms for solving the inverse problems of model parameters estimation require a numerical solution of the corresponding direct problem.Moreover, if an inverse problem reduces to the optimization problem, the corresponding direct problem should be solved many times.However, finding a numerical solution of the direct problem for the TFDPL model is computationally expensive due to memory effects.
The algorithm utilizes a small amount of experimental data.The temperature measurement should be performed only in two points: on the surface of the body and at any interior point.Hence, only two thermal sensors are needed.This reduces the cost of experimental research.
The algorithm uses explicit representations for three model parameters -the thermal diffusivity, a, the fractional analogue of thermal relaxation time, τ q , and the fractional analogue of thermal retardation, τ T .All integrals in these representations can be efficiently evaluated by quadratures.The direct Laplace transform can be made numerically by using the Gauss-Laguerre quadrature (see, e.g., [49]).Thus, the algorithm is simple and has a small time complexity.
There are only two internal parameters in the algorithm: p min and p max .These parameters can be considered as regularization parameters and allows one to obtain a stable solution of inverse problem.
However, the present version of the algorithm is limited to a one-dimensional linear TFDPL model with constant parameters for a half-space heat transfer medium.
Numerical calculations show (see results in Section 5) that for a given fractional order α the algorithm permits one to estimate other parameters of the TFDPL model with a good accuracy if the parameters p min and p max are chosen correctly.The value of p max is controlled by the determinant of the matrix A from the system (38).For stability, it is recommended that det(A) = O(10 k ) with k ≥ −3.The value of p min is highly depended on α: the decreasing of α leads to the decreasing of p min , so that p min ≈ p 1/α 0 , where p 0 = p min | α=1 < 1.In practice, the value of p min can be obtained iteratively from the condition ψ(p min ) ≈ a −1 , where a is the thermal diffusivity and the function ψ(p) is defined in (36).In this case, the estimation errors and errors of input data have the same order of magnitude.
The proposed algorithm permits one to estimate all four parameters of the TFDPL model only if the relative error of input data is less than 1%.Despite the fact that thermal sensors can measure the temperature with much less errors, this is a strict limitation of the algorithm.The cause of this limitation is that Equation ( 35) is nonlinear with respect to fractional order α and linear with respect to other parameters of the TFDPL model.As a result, even a small perturbation of α leads to a high perturbation of a, τ t , and τ q .Note that τ q and b = aτ T are multipliers in the terms with p α , whereas a is not one.For this reason, the thermal diffusivity a is more sensitive to a perturbation of α.Since τ T = b/a, its error depends on the error of a.The example given in Section 5 demonstrates these trends.The fractional order α governs the asymptotic behaviour of the temperature field for large times.Therefore, it can be potentially estimated from this asymptotic.However, how we can incorporate this idea into the TIC algorithm is still an open question.
Finally, we will pay attention to another important issue.As mentioned in [27], confidence intervals for estimated parameters should be established during a parameter estimation process.An example of the successful calculation of such intervals can be found in [58].However, the proposed algorithm does not allow one to compute these intervals in the classical manner.This is because we work in the Laplace space and deal with the continuous function ψ(p).Thus, it is a significant task to develop an approach for evaluating the confidence intervals of the parameters that were estimated using the proposed algorithm.

Conclusions
In this paper, we propose a new semi-explicit algorithm for parameters estimation in a time-fractional dual-phase-lag heat conduction model.The algorithm was developed by using the time integral characteristics method and allows one to identify the model parameters using the temperature measurements at only one interior point.The algorithm has a low computational cost because it does not require a numerical solution of the corresponding direct heat conduction problem.This makes it different from most other algorithms of parameters estimation.We provide the recommendations for choosing numerical values of the internal parameters of the algorithm.A numerical example illustrates its strengths and weaknesses.The algorithm was obtained under the assumption of half-space heat conduction medium.In particular, this assumption is fulfilled during laser heating.However, overcoming this assumption is a possible direction of future work.The algorithm showed a low accuracy for estimating the order of fractional differentiation if temperature measurement errors are high enough.Solving this problem is another important direction for further research.Finally, developing an approach for calculating the confidence intervals of the estimated parameters is also a challenging problem.

Figure 1 .
Figure 1.Graphs of the function ψ(p) defined by (42) for different values of α.

Figure 5 .
Figure 5. Graphs of det(A) for a fixed p min and various p max .

Figure 6 .
Figure 6.Graphs of det(A) for a fixed p max and various p min .

Table 1 .
Comparison of the restored parameters for different approximations of ψ(p) with α = 0.75.

Table 2 .
Comparison of the restored parameters including α for different approximations of ψ(p).

Table 3 .
Comparison of the restored parameters for different approximations of ψ(p) with α = 0.5.

Table 4 .
Comparison of the restored parameters including α for different approximations of ψ(p).

Table 5 .
Comparison of the restored parameters for different approximations of ψ(p) with α = 0.25.Table6contains the results of parameters estimation for unknown α.It can be seen from the table that, in all cases, we have a high level of error, especially for τ q .

Table 6 .
Comparison of the restored parameters including α for different approximations of ψ(p).

Table 7 .
Comparison of the restored parameters for different approximations of ψ(p) with α = 0.25.

Table 8 .
Comparison of the restored parameters including α for different approximations of ψ(p).