Long-memory recursive prediction error method for identification of continuous-time fractional models

This paper deals with recursive continuous-time system identification using fractional-order models. Long-memory recursive prediction error method is proposed for recursive estimation of all parameters of fractional-order models. When differentiation orders are assumed known, least squares and prediction error methods, being direct extensions to fractional-order models of the classic methods used for integer-order models, are compared to our new method, the long-memory recursive prediction error method. Given the long-memory property of fractional models, Monte Carlo simulations prove the efficiency of our proposed algorithm. Then, when the differentiation orders are unknown, two-stage algorithms are necessary for both parameter and differentiation-order estimation. The performances of the new proposed recursive algorithm are studied through Monte Carlo simulations. Finally, the proposed algorithm is validated on a biological example where heat transfers in lungs are modeled by using thermal two-port network formalism with fractional models.


Introduction
Fractional calculus was explored in the twentieth century as a way to model diffusion phenomena for different geometries [32,33]. It has been shown that diffusion phenomena can be modeled by using fractionalorder transfer functions with differentiation orders that are multiples of 0.5. In electro-chemistry, charge diffusion in acid batteries is governed by Randles model, which involve a half-order fractional integrator [38]. Studies have also shown a half-order fractional integrator model for semi-infinite thermal systems [1]. These physical phenomena of diffusion have led to the theoretical conception of constant-phase elements [16,30] that may be used as building blocks for circuit models. Constant phase elements may be used to model intestine tissue behavior [9], cardiac tissue [21], porous films [4] and lung impedance behavior [14,44].
New challenges appear on system identification: [46] provides new paradigms and challenges in system identification such as broader types of uncertainties, networked systems or even data explosion; [35] proposes kernel methods; [20] gives new kernel-based regularization methods; [2] proposes to estimate time delay with sampled limit cycle in frequency domain, etc. There is a vast literature concerning the analysis and identification of discrete-time models [19]. However, note that such type of models may lead to parameters that lack physical meaning. The sampling zero problem [48] related to discretization should not be neglected either. As a model with physical meaning is sought, a continuous-time identification is better suited. The main challenge with adapting system identification for continuous-time models from sampled data relies on the differentiation approximations for the information matrix [10]. When dealing with fractional models, this information matrix will contain approximations of fractional-order differentiations from experimental data.
Since the 1990s, substantial contributions have been made for fractional-order system identification by estimating the coefficients of a fractional-order transfer function. In continuous-time identification, methods relying on state variable filters, least squares and instrumental variables were introduced in [3,31,42]. These techniques have been applied for identifying thermal diffusion on an aluminum rod [24,42], charge diffusion lithium-ion batteries [8], muscle relaxation [40], heat transfers in biological system [44].
Most of the fractional system identification methods work offline, which means that the whole data of an experiment are used to provide parameter estimation. In some cases, online parameter estimation is needed such as in open-heart surgery where lung thermal online system identification will prevent pulmonary cell dying. The most famous online identification method is the recursive least squares. Recursive prediction error methods are also available in the literature [19]. Padilla [34] provides a recursive version of LSSVF and RIVC for integer-order systems. By starting from Le Lay's method and online estimation concepts, a method to estimate coefficients of a fractionalorder transfer function was proposed in [6]. Fractionalorder models may also require the estimation of one or several differentiation orders, which has already been studied in the offline cases by means of gradient descent for single-input single-output fractional systems [42] or even in the multiple-input single-output fractional systems [26,45]. An extension of Djouambi's identification method has also been proposed to identify the fractional orders by using a scheme similar to gradient descent in [13].
The aim of this paper is to first provide efficient recursive all parameter estimation methods, namely both the coefficients and the differentiation orders. First of all, when the differentiation orders are assumed known, a new method called long-memory recursive prediction error method is proposed and compared to existing methods through Monte Carlo simulation analysis. A specific aspect of fractional differential equation modeling is the determination of the differentiation orders. In system identification of classic rational models, as the model order remains unchanged, only the coefficients are estimated. When estimating both coefficients and differentiation orders of fractional models, this task is not trivial: Indeed, in an online recursive algorithm, the model order changes at each iteration. Therefore, the long-memory prediction error method (LMRPEM) is extended for estimating both parameters and differentiation orders. Monte Carlo simulation analysis is again proposed to analyze the statistical properties of our new method. Finally, the new LMR-PEM algorithm is applied for identifying heat transfers on a human bronchus.
Real experiments for modeling heat transfers in sheep lung were previously proposed in [44] for system identification by using black box models through Havriliak-Negami functions. It is now proposed to use a more realistic model, which is closer to the physics that governs heat transfers in human lungs. By using thermal two-port network formalism [22], which is none other than the solving of the heat equation on a unidirectional medium, and fractional-order impedance approximations [7], heat transfers in lungs of a human main bronchus can be modeled as an equivalent T network (see Fig. 5). All analytical developments for heat transfer modeling are detailed in [7].
The paper is organized as follows: Sect. 2 describes the problem formulation for recursive fractional system identification, Sect. 3 presents the recursive estimation methods for coefficient estimation of fractional systems, and then, Sect. 4 introduces a recursive estimation method for all-parameter estimation by combining both coefficient and fractional differentiation-order estimation. Section 5 presents an application to thermal transfer in human lungs, and conclusions are drawn in Sect. 6.

Mathematical background
Many definitions exist for fractional differentiations [11,15,39]. Grünwald [12] and Letnikov [17] independently developed a definition for an arbitrary-order differentiation in terms of a convergent series: where . stands for the floor operator, p is the differential operator p = d dt , and α j is Newton's binomial coefficient generalized to real numbers: where is Euler's Gamma function. The slow convergence of this binomial justifies why fractional-order operators are well suited for long-memory behavior modeling, as the whole past of the function is required for its definition. By replacing h by the sampling time T s , an implementable approximation of the fractional differentiation is obtained: Note that the error terms are proportional to the sampling time [36,Section 7.4]. Consequently, the smaller the sampling time, the lesser the approximation errors. 1 If one considers a case with null initial conditions, the Laplace transform of the fractional differentiation may lead to a simple expression [36,Section 2.8.4]: For a fractional single-input single-output (SISO) model, a fractional differential equation relates the input u(t) to the output y(t): where (a i , b j ) are real numbers and α 1 < α 2 < . . . < α N and β 0 < β 1 < . . . < β M are non-integer positive numbers.
The Laplace transform of (5) gives the following fractional-order transfer function model: Definition 1 A SISO system G(s) (6) is commensurate of order ν, if all differentiation orders are integer multiples of ν: where all j = β j ν and i = α i ν are integer numbers.
Stability of fractional-order systems has been analyzed in different contexts. The most well-known stability criterion was established by Matignon [25] and allows to check the stability of a commensurate-order system through the location of its s ν -poles. The original theorem was established for commensurate orders ν, for 0 < ν < 1, but was extended for orders between 1 and 2 in [29]. Orders with commensurate orders beyond 2 are unstable (see [23]). Further extensions have been developed in order to check stability of noncommensurate systems [37].

Theorem 1 Let S be a commensurate transfer function and ν its commensurate order. G(s) = Q ν (s) P ν (s) is BIBOstable if and only if:
and, for every pole s k (P ν (s k ) = 0):

Problem formulation
Let be an input u(t) and its noise-free output y(t) related by equation (6), where A(s) and B(s) are assumed coprime and the system stable. Data are collected at regular samples from t = 0 to the current time t = kT s with the sampling time T s . Furthermore, the number of data samples is assumed to be large enough so that the convergence of the estimated parameters to the true ones is guaranteed. Moreover, the quasistationary input signal {u(t), 0 ≤ t ≤ kT s } applied to the system is persistently exciting and gives rise to an output signal {y(t), 0 ≤ t ≤ kT s }. The output data are corrupted by an additive white measurement noise ξ(t) due to experiment and sensor imperfections. The noise is assumed normally distributed with a zero mean and variance, σ 2 , considered at discrete instants. The complete system equation can be written as: where y * (t) is the noisy sampled value of the unobserved noise free output y(t). When using model (6), the parameter vector, θ , is defined as where ρ is the vector of N + M + 1 coefficients: and μ is the vector of N + M +1 differentiation orders: The total number of parameters to be estimated is then 2×(N + M +1). If N or M are too high, the complexity may be too high and nonlinear optimization algorithms may fail to converge to a global minima.
In order to reduce the number of parameters, a commensurate-order model (7) can be used, where the differentiation order vector μ is reduced to a single parameter ν. The parameter vector θ is then reduced to N + M + 1 coefficients as in (12) added up of the commensurate order μ = ν.

Coefficient estimation of fractional-order systems
In this section, all differentiation orders are assumed known; therefore, only the coefficients are estimated. The parameter vectorθ boils down to: 3.1 Recursive least squares with state variable filter (RLSSVF) The estimated output can be expressed in a linear form as: where the regression vector is as follows A suitable error function is defined as: where the parameter vectorθ minimizes the quadratic error function: which can be put under the form of an optimization problem: It is well known that the least squares are biased in the presence of noise, and therefore, state variable filters (SVFs) can be used to obtain a filtered input u f (t) and output y f (t). This filtering allows to reduce the noise influence on the estimation and is normally a lowpass filter tuned to emphasize a particular bandwidth. A typical choice is the low-pass Poisson filter: The filter should be designed by respecting the following [3]: where · is the ceil operator. It has been shown in [5] that the cutoff frequency ω f may have an influence on the convergence rate of the algorithm, as well as sampling time and input. This frequency is typically chosen to be somewhat larger to the system bandwidth [34]. This means that a priori estimation or tuning of the system bandwidth may be required to get better results.
The error function (t) is now replaced by its filtered version f (t): where the filtered regression vector φ * . (23) and the filtered input u f (t) and output y * f (t) are defined as: Inspired by the recursive least squares [19, chap. 11], the RLSSVF algorithm is obtained for fractional-order systems: where γ is a forgetting factor that can give more or less weight to the most recent past of signals and F is the covariance matrix. Remark: the forgetting factor γ may actually be a function of time and be adaptive as identification is running [41]. In this study, it will be taken as a constant and equal to unity.

Recursive prediction error method (RPEM)
By following the prediction error method, a suitable error function is given by the output error: where the estimated outputŷ(t) is computed as: The quadratic error function defined in (18) is minimized. By taking the gradient of the criterion: with its error sensitivity is computed by: where and Inspired by the recursive least squares [19, chap. 11], the RPEM algorithm is obtained for fractional-order systems: where γ ρ is a refining gain analogous to the step in a gradient descent. This forgetting factor γ ρ may actually be adapted during the identification or be a function of time. However, the determination of an adequate γ ρ may prove to be non-trivial, as stated in [18]. For the present work, it will be taken as a constant.

Contribution with long-memory recursive prediction error method (LMRPEM)
Fractional-order systems have a natural long-memory behavior as the whole past of a signal differentiation is needed, and therefore, an extended error with long memory is more suited: This error signal˜ will include errors at all instants from t = 0 to the current time t = kT s and will be increased by one data-point per iteration. Consequently, the resulting gradientψ ρ (kT s ) is now a matrix: where k indicates the present iteration. By introducing the extended measured outputỸ * (kT s ) and the extended estimated outputỸ (kT s ) as: the long-memory recursive prediction error method (LMRPEM) is proposed for fractional-order systems: where again, γ ρ is a refining gain analogous to the step in a gradient descent. As previously stated in Sect. 3.2, γ ρ will be taken as a constant.

True system description
Let us consider the following true system model: with K = 1, ζ = −0.5, ω 0 = 0.5rad/s and ν = 0.5. Such a fractional model is called a second kind fractional-order system [23] and has a resonant peak around ω r = 0.45rad/s (see Fig. 1). In fact, this is a typical kind of model structure that well represents diffusive phenomenon such as thermal diffusion [43]. Moreover, such a model structure is similar to the one that models diffusive phenomenon in lungs (see the human bronchus heat transfer application in Sect. 5).
The model structure is chosen such as the true one (39) and is expressed by: with differentiation orders assumed known and fixed just like the true system (39), where ν = 0.5. The aim is to estimate the system coefficients and to see the statistical properties of the proposed LMRPEM algorithm through Monte Carlo analysis, with N exp = 100 runs and different signal-to-noise (SNR) levels.
The input signal is a pseudo-random binary sequence oscillating between − 5 and 5 and containing 1000 values with a sampling time T s = 0.01s. All simulations  Figure 2 shows the input used for simulations, as well as the noise-free output and the noisy output with a SNR = 15 dB. All simulations have been carried out with the CRONE toolbox developed in MATLAB , which is dedicated to fractional calculus, fractional system simulation and system identification with fractional models. 2

Coefficient estimation with known differentiation orders
From Fig. 1, the true system is resonant in the frequency range [1, 10]rad/s. Therefore, the Poisson filter (20) is set with a cutoff frequency of ω f = 10rad/s and as max(β M , α N ) = 1, N f = 2. Note that in a realistic scenario such information may not be available and an empirical estimation of the cutoff frequency ω f and the order N f will be required.
The Monte Carlo simulation results for 100 runs are summarized in Table 1.
The RLSSVF method lacks robustness with respect to measurement noise. The lesser the SNR is, the greater the estimation bias gets and the greater the standard deviations become. Moreover, the RLSSVF is not efficient as the parameters are not correctly estimated for any noise level. An improvement can be obtained by adjusting the cutoff frequency of the filter; however, it should be noted that a frequency ω f extremely close to the system bandwidth would deteriorate the estimation in a practical scenario.
The RPEM method (γ ρ was set to 0.01) is more consistent than the RLSSVF method as the bias is less present. However, the lesser the SNR, the greater the bias. The true value of the parameters is always included in the range determined by the standard deviation. On the other hand, a rough factor of 3 can be observed between the standard deviations of SNR = 20 dB and 10 dB.
The LMRPEM method (γ ρ was set to 0.01) has the most consistent estimation as the estimation has no bias and the estimation standard deviations (and consequently the estimation variances) are the lowest whatever the noise level. This method has significantly reduced the parameter deviations.

Commensurate order with unknown differentiation orders
In practice, the differentiation orders are not always known a priori. The commensurate-order influence is evaluated by computing the cost function, defined as the 2 -norm of the normalized output error: where the long-memory error˜ (N T s ,θ) is defined as in (34): from relations (36) and (37).  Figure 3 shows the criterion J versus the commensurate order obtained by using LMRPEM method for different values of the commensurate order in order to approximate the true system with a noise-to-signal ratio N S R = −20 dB. As expected, the criterion minimum is obtained at ν = 0.5, which is the true commensurate order. Indeed, when the model is in the same model class as the true system, the minimum of the cost function is found at ν = 0.5 and equals − 20 dB which exactly corresponds to the NSR. Other values of ν lead to non-negligible modeling errors as the coefficient estimation would not converge to the true parameters. For example, at ν = 1, the cost function is around − 10.01 dB, which corresponds to a modeling error around 10 dB.
This simulation result motivates to estimate the fractional differentiation orders, as they may considerably influence the estimation.

Differentiation-order estimation and all-parameter system identification
In this section, the differentiation orders are assumed unknown, as it is often the case in practice. Thus, it is helpful to consider differentiation-order estimation along with coefficient estimation. As the differentiation orders are not expressed in a linear way in (6) or (7), it is necessary to recursively estimate these parameters with a suitable error function. As the LMRPEM method has given the best estimation results in the pre- vious section for coefficient estimation, it is proposed to use this method for the differentiation-order estimation. Then, all parameters are estimated for both, the coefficients and the differentiation-order estimation, through a two-stage algorithm.

Differentiation-order estimation
In a general context, a vector of differentiation orders could be determined and when dealing with a commensurate model, a single differentiation order is estimated. By following the prediction error method, a suitable error function is given by the output error: and differentiation orders may be recursively adjusted by adding a correction term to the previous iteration: where Recall that the LMRPEM method has given the best results for estimating the coefficients for fractionalorder systems. Consequently, an extended error is used to take into account the long-memory behavior of the system: (46) which turns the error sensitivity into a matrix: For commensurate-order estimation, the gradient boils down to: The error sensitivity function is defined by: where the commensurate-order output sensitivity ∂ŷ ∂ν is expressed by: As time simulation of such a function is non-trivial, the gradient will be numerically computed and the LMR-PEM will be used for differentiation-order estimation: 4.2 All-parameter system identification with a hybrid two-stage algorithm The best previous method for coefficient estimation, namely the LMRPEM, is combined with differentiation-order estimation, through another LMRPEM. This two-stage algorithm, called LMRPEM-2, is used for all-parameter estimation of a fractional-order model.

38) (b) Differentiation order estimation
Assuming the coefficient ρ(k) known, compute ν(k) from equation (51) As stated by Young in [47], the tuning of the forgetting factor γ in recursive PEM methods is non-trivial and influences the convergence of the parameters. Consequently, the LMRPEM-2 algorithm 1 requires the tuning of two forgetting factors (one for the coefficients and one for the differentiation orders).

Simulation results for recursive all-parameter estimation
The same system (39) as detailed in Sect. 3.4.1 is used with different noise levels: ∞ without noise, 20 dB, and 15 dB. The tuning of the γ gains was performed empirically and is given in Table 2.  The parameter evolution for an SNR = 20 dB is plotted in Fig. 4 for Algorithm 1. On the commensurateorder estimation, hybrid LMRPEM-2 Algorithm 1 exhibits some oscillations at the beginning. However, the hybrid LMRPEM-2 algorithm tends to the true commensurate order ν = 0.5.
Monte Carlo simulation analysis is provided in Table 3. Three levels of noise are considered, and for each level, 100 runs have been carried out. Note that all runs have converged to stable models. In noiseless conditions, the hybrid LMRPEM-2 Algorithm 1 has converged to the true commensurate order, and therefore, the coefficient estimations are pretty accurate in mean. As the noise level increases, hybrid LMRPEM-2 Algorithm 1 presents consistent estimates, with low variance and without bias.
The present work shows parameter estimation within a finite data length. In practice, realistic applications may require much longer datasets. One main issue regarding longer datasets lies on the computation time. Grünwald-Letnikov's differentiation is longer to calculate at each iteration when the dataset gets long enough, as the whole past of the signal is used for its computation. If the dataset gets even longer, the computation time will definitely be affected. Therefore, an extremely long dataset would require an extremely long calculation time by the end of the simulation, which would be impossible to implement in a real-time scenario. By allowing some approximation errors, the past of the function used may be limited to a memory length L, as stated by the short memory principle [36].

Application: human bronchus heat transfer system identification
Medical and surgery applications need online biological parameter estimation: Time can be highly critical, and more especially, accurate real-time estimation is required. Real experiments for modeling heat transfers in sheep lung were previously proposed in [44] for system identification by using black box model through Havriliak-Negami functions. It is now proposed to use a more realistic model for heat transfers in human lungs. By using thermal two-port network formalism [22], which is none other than the solving of the heat equation on a unidirectional medium, and fractional-order impedance approximations [7], heat transfer in lungs of a human main bronchus can be modeled as an equivalent T network (see Fig. 5). All model developments are detailed in [7]. The bronchus is modeled as an air cylinder with an intermediate length of L = 0.0236 m [28] and radius r ≈ 1 mm [27]. Impedance approximations are given by: and For an insulated output case (φ out = 0), the transfer function relating the output temperature to the heat flux The input signal is a pseudo-random binary sequence oscillating between − 1 and 1 and containing 2000 values with a sampling time T s = 1s. The initialization is θ T = [2, 5, 5, 10] and ν = 0.35. The noise was set to SNR = 20 dB. The input/output data are shown in Fig. 6.
This model is set as the true system (54): Hybrid LMRPEM-2 Algorithm 1 is used to estimate the parameters as it is the most efficient algorithm for all-parameter estimation (see Sect. 4.3). γ ρ and γ ν are the same as given in Table 2. The recursive coefficient and differentiation-order estimations are shown in Fig.  7. The commensurate-order convergence is slow but  Fig. 7 Parameter estimation for bronchus application very accurate toward the system true commensurate order ν = 0.5. The coefficient convergence exhibits fairly slow convergence rate, especially for a 2 ; however, it tends to its true value toward the end of the simulation. Therefore, when the dataset gets long enough, the convergence to the true parameters is improved. Finally, at the last iteration, the estimated model is expressed as the following: G(s) = 2.519s 0.502 + 1.013 12.66s 1.005 + 2.601s 0.502 + 1 , an estimated model that has well converged to the true parameters such as predicted by the hybrid LMRPEM-2 Algorithm 1.
Note that the coefficients have correctly converged to their true values only after the convergence of the commensurate order. Consequently, well estimating the differentiation orders is necessary so that the coefficients well converge, and therefore, physical parameters can be extracted.

Conclusions and perspectives
A new method called LMRPEM-2 is proposed for allparameter recursive system identification of continuous-time fractional-order models. Least squares (RLSSVF) and prediction error method (RPEM) were compared to a new variant of the RPEM, called LMRPEM. The latter method takes into account the long-memory behavior of fractional-order systems. RLSSVF method showed poor performances, as it provides high bias as noise increases and implies a priori estimation of the system bandwidth. Other methods provided more accurate coefficient estimation as they present no bias. The statistical properties of the LMRPEM have proven to produce the lesser parameter variance as compared to the other methods.
Then, for all-parameter estimation, especially for differentiation order estimation, a new hybrid algorithm is introduced: the hybrid LMRPEM-2 for both coefficient and differentiation-order estimation. The LMRPEM-2 algorithm proved to be consistent: Parameter estimations are without bias with very low variance. Also, the longer the dataset, the better the convergence to the true parameters as fractional systems have long-memory behavior.
Finally, a biological example was proposed: Heat transfers in lungs are modeled on a human main bronchus by using thermal two-port network formalism with fractional models. Applying the hybrid LMRPEM-2 method has provided consistent estimates without bias and with very low variance.
Research perspectives may include further studies regarding a method to tune the forgetting factor γ in a methodical way. The issue regarding computation time with long datasets should also be studied by taking into account the short memory principle as well as the frequency content of the analyzed signals. Real experiments on sheep will be carried out at the Bordeaux University Hospital and IHU Liryc Institute for online system identification of lung heat transfers.
Author contributions J.-F. Duhé has developed the recursive algorithms with S. Victor. He also made the simulations. S. Victor, P. Melchior, Y. Abdelmounen and F. Roubertie have contributed to the biological and physiological results.

Funding Not applicable.
Data availability Enquiries about data availability should be directed to the authors.