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

Stephane Victor, Jean Francois Duhe, Pierre Melchior, Youssef Abdelmounen, François Roubertie

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Stephane Victor, Jean Francois Duhe, Pierre Melchior, Youssef Abdelmounen, François Roubertie. Long-memory recursive prediction error method for identification of continuous-time fractional models. Nonlinear Dynamics, 2022, 110 (1), pp.635-648. 10.1007/s11071-022-07628-8. hal-03865862

HAL Id: hal-03865862

https://hal.science/hal-03865862

Submitted on 22 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Long Memory Recursive Prediction Error Method for Identification of Continuous-time Fractional Models

**stephane victor** (✉ stephane.victor@ims-bordeaux.fr)

IMS: Laboratoire de l'integration du materiau au systeme https://orcid.org/0000-0002-0575-0383

**Jean-François Duhé**

IMS: Laboratoire de l'integration du materiau au systeme

**Pierre Melchior**

IMS: Laboratoire de l'integration du materiau au systeme

**Youssef Abdelmoumen**

Clinique Saint Augustin

**François Roubertie**

INSERM U1045

## Research Article

**Keywords:** continuous-time models, fractional calculus, fractional order model, system identification, recursive identification, real-time system identification, instrumental variable, prediction error method, least squares, long memory prediction error method

**Posted Date:** January 24th, 2022

**DOI:** https://doi.org/10.21203/rs.3.rs-1272889/v1

**License:** © This work is licensed under a Creative Commons Attribution 4.0 International License.

Read Full License
# Long memory recursive prediction error method for identification of continuous-time fractional models

Stéphane Victor · Jean-François Duhé · Pierre Melchior · Youssef Abdelmounen · François Roubertie

Received: date / Accepted: date

**Abstract** This paper deals with recursive continuous-time system identification using fractional differentiation 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 of the classic methods used for integer order models, are compared to our new method, thus proving the efficiency of our algorithms. 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 is studied through Monte-Carlo simulations. Finally, the proposed algorithm is validated on a biological example where heat transfer in lungs is modeled by using thermal two-port network formalism with fractional models.

**Keywords** continuous-time models · fractional calculus · fractional order model · system identification · recursive identification · real-time system identification · instrumental variable · prediction error method · least squares · long memory prediction error method

## 1 Introduction

Fractional calculus was explored in the twentieth century as a way to model diffusion phenomena for different geometries [31],[32]. It has been shown that diffusion phenomena can be modelled by using fractional order 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 [37]. Studies have also shown a half-order fractional integrator model for semi infinite thermal systems [1]. This physical phenomena have led to the theoretical conception of constant-phase element [16,29], that may be used as building blocks for circuit models. Constant-phase element may be used to model intestine tissue behavior [9], cardiac tissue [21], porous films [4] and lung behavior [14,42].

New challenges appear on system identification: [43] provides new paradigms and challenges in system identification such as broader types of uncertainties, networked systems or even data explosion; [34] proposes kernel methods; [20] gives new kernel-based regularization methods; [2] proposes to estimated 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 [45] 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 in the derivative approximations for the information matrix [10]. When dealing with fractional models, this information

<footnote>S. Victor, J.-F. Duhé, and P. Melchior
Univ. Bordeaux, CNRS, IMS UMR 5218, Bordeaux
INP/enseirb-matmeca, 351 cours de la Libération, 33405
Talence Cedex France Tel.: +33 5 40 00 36 27 E-mail: stephane.victor@ims-bordeaux.fr
Y. Abdelmounen and F. Roubertie
IHU Liryc, Electrophysiology and Heart Modeling Institute, Fondation Bordeaux Université 33000 Bordeaux Cedex
France
INSERM U1045, Centre de recherche Cardio-Thoracique de
Bordeaux, 33000 Bordeaux Cedex France
Y. Abdelmounen
Clinique Saint Augustin</footnote>
matrix will contain approximations of fractional order derivatives 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, 41, 30]. These techniques have been applied for identifying thermal diffusion on an aluminium rod [24, 41], charge diffusion lithium-ion battery models [8], muscle relaxation [39], heat transfers in biological system [42].

Most of the fractional system identification methods work offline, which means that the whole data of an experiment is used to provide parameter estimation. In some cases, online parameter estimation is needed in open-heart surgery where lung thermal modeling will prevent pulmonary cell dying. The most famous online identification method is the recursive least squares. Recursive prediction-error methods are also available in literature [19]. Padilla [33] provides a recursive versions 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 fractional order transfer function was proposed in [6]. Fractional order 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 [41]. An extension of Djouambi's identification method has also been proposed to identify the order 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 analysis the statistical properties of our new method. Finally, the new LMRPEM algorithm is applied for identifying heat transfers on a human bronchus.

Real experiments for modeling heat transfers in sheep lung was previously proposed in [42] 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 transfers in lungs of a human main bronchus can be modeled as an equivalent T network (see figure 5). All analytical developments for heat transfer modeling are detailed in [7].

The paper is organized as follows: Section 2 describes the problem formulation for recursive fractional system identification, Section 3 presents the recursive estimation methods for parameter estimation of fractional systems, then Section 4 introduces 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 Section 6.

## 2 Problem formulation for recursive fractional system identification

### 2.1 Mathematical background

Many definitions exist for fractional derivatives [15, 11, 38].

Grünwald [12] and Letnikov [17] independently developed a definition for an arbitrary order derivative in terms of a convergent series:

$$p^{\alpha} f(t) = \lim_{h \to 0} \frac{1}{h^{\alpha}} \sum_{j=0}^{\lfloor \frac{t}{h} \rfloor} (-1)^j \binom{\alpha}{j} f(t - jh) \quad (1)$$

where $\lfloor \cdot \rfloor$ stands for the floor operator, $p$ is the differentiation operator $p = \frac{d}{dt}$, and $\binom{\alpha}{j}$ is Newton's binomial coefficient generalized to real numbers:

$$\binom{\alpha}{j} = \frac{\Gamma(\alpha + 1)}{\Gamma(j + 1)\Gamma(\alpha - j + 1)} = \frac{\alpha(\alpha - 1)\cdots(\alpha - j + 1)}{j!} \quad (2)$$

where $\Gamma$ 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$, one gets an implementable approximation for the fractional derivative.
tive:

$$p^{\alpha} f(t) = \frac{1}{T_s^{\alpha}} \sum_{j=0}^{\lfloor \frac{t}{T_s} \rfloor} (-1)^j \binom{\alpha}{j} f(t - jT_s) + O(T_s). \quad (3)$$

Note that the error terms are proportional to the sampling time [35, Section 7.4]. Consequently, the smaller the sampling time, the lesser the approximation errors<fnref n="1" />.

If one considers a case with null initial conditions, the Laplace transform of the fractional derivative may lead to a simple expression [35, Section 2.8.4]:

$$\mathcal{L}\{p^{\alpha} f(t)\} = s^{\alpha} F(s). \quad (4)$$

For a fractional SISO model, a fractional differential equation relates the input $u(t)$ to the output $y(t)$:

$$y(t) + a_1 p^{\alpha_1} y(t) + \dots + a_N p^{\alpha_N} y(t) = b_0 p^{\beta_0} u(t) + b_1 p^{\beta_1} u(t) + \dots + b_M p^{\beta_M} u(t) \quad (5)$$

where $(a_i, b_j)$ are real numbers and $\alpha_1 < \alpha_2 < \dots < \alpha_N$ and $\beta_0 < \beta_1 < \dots < \beta_M$ are non-integer positive numbers.

The Laplace transform of (5) gives the following fractional order transfer function model:

$$G(s) = \frac{B(s)}{A(s)} = \frac{\sum_{i=0}^{M} b_i s^{\beta_i}}{1 + \sum_{j=1}^{N} a_j s^{\alpha_j}}. \quad (6)$$

**Definition 1** A SISO system $G(s)$ (6) is commensurate of order $\nu$, if all differentiation orders are integer multiples of $\nu$:

$$G(s) = \frac{B(s)}{A(s)} = \frac{\sum_{i=0}^{M} b_i s^{i\nu}}{1 + \sum_{j=1}^{N} a_j s^{j\nu}}, \quad (7)$$

where all $j = \frac{\beta_j}{\nu}$ and $i = \frac{\alpha_i}{\nu}$ 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^\nu$-poles. The original theorem was established for commensurate orders $0 < \nu < 1$, but this was extended for orders between 1 and 2 [28]. Orders with commensurate orders beyond 2 can be proven to be unstable [23]. Further extensions have been developed in order to check stability of non-commensurate systems [36].

**Theorem 1 (Matignon's stability Theorem)** Let $S$ be a commensurate transfer function and $\nu$ its commensurate order. $G(s) = \frac{Q_{\nu}(s)}{P_{\nu}(s)}$ is BIBO-stable if and only if:

$$0 < \nu < 2 \quad (8)$$

and, for every pole $s_k$ ($P_\nu(s_k) = 0$):

$$|\arg(s_k)| > \nu \frac{\pi}{2}. \quad (9)$$

## 2.2 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 current time $t = T_f$ with the sampling time $T_s$. Furthermore, the number of data sample is assumed to be large enough so that the convergence of the estimated parameters to the true ones is guaranteed. Moreover, the quasi-stationary input signal $\{u(t), 0 \le t \le T_f\}$ applied to the system is persistently exciting, and gives rise to an output signal $\{y(t), 0 \le t \le T_f\}$. The output data is corrupted by an additive white measurement noise $\xi(t)$ due to experiment and sensor imperfections, normally distributed with a zero mean and $\sigma^2$ variance, considered at discrete instants. The complete equation can be written as:

$$S: \begin{cases} y(t) = G(\rho)u(t) \\ y^*(t_k) = y(t_k) + \xi(t_k), \end{cases} \quad (10)$$

where $y^*(t_k)$ is the noisy sampled value of the unobserved noise free output $y(t)$.

When using model (6), the parameter vector, $\theta$, is defined as

$$\theta = [\rho \quad \mu]^T \quad (11)$$

where $\rho$ is the vector of $N + M + 1$ coefficients:

$$\rho = [b_0, b_1, \dots b_N, a_1, a_2, \dots a_M]^T, \quad (12)$$

and $\mu$ is the vector of $N + M + 1$ differentiation orders:

$$\mu = [\beta_0, \beta_1, \dots \beta_M, \alpha_1, \alpha_2, \dots \alpha_N]^T. \quad (13)$$

The total number of parameters to be estimated is then $2 \times (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 $\mu$ is reduced to a single parameter $\nu$. The parameter vector $\theta$ is then reduced to $N + M + 1$ coefficients as in (12) added up of the commensurate order $\mu = \nu$.

<footnote n="1">¹ If the sampling time is too small, numerical problems may occur such as stability in digital implementation. In this case, suitable discrete rational approximations could be used.</footnote>
### 3 Coefficient estimation of fractional order system

In this section, all differentiation orders are assumed known, therefore, only the coefficients are estimated. The parameter vector $\hat{\theta}$ boils down to:

$$ \hat{\theta} = \hat{\rho} = [\hat{b}_0, \hat{b}_1, \dots \hat{b}_M, \hat{a}_1, \hat{a}_2, \dots \hat{a}_N]^T. \quad (14) $$

#### 3.1 Recursive Least-Squares with State Variable Filter (RLSSVF)

The estimated output can be expressed in a linear form as:

$$ \hat{y}(t) = \phi^*(t)^T \hat{\theta} \quad (15) $$

where the regression vector is as follows

$$ \phi^*(t)^T = \begin{bmatrix} p^{\beta_0} u(t), \dots, p^{\beta_M} u(t), \\ -p^{\alpha_1} y^*(t), \dots, -p^{\alpha_N} y^*(t) \end{bmatrix}. \quad (16) $$

A suitable error function is defined as:

$$ \epsilon(t) = y^*(t) - \phi^*(t)^T \hat{\theta}, \quad (17) $$

where the parameter vector $\hat{\theta}$ minimizes the quadratic error function:

$$ J(\hat{\theta}) = \frac{1}{2} \sum_{k=1}^{t} \epsilon(k, \hat{\theta})^2, \quad (18) $$

which can be put under the form of an optimization problem:

$$ \hat{\theta}(t) = \arg \min \sum_{k=1}^{t} \epsilon(k)^2. \quad (19) $$

It is well-known that the least-squares are biased in presence of noise, therefore, state-variable filters 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 low-pass filter tuned to emphasize a particular bandwidth. A typical choice is the low-pass Poisson filter:

$$ F(p) = \frac{1}{\left(1 + \frac{p}{\omega_f}\right)^{N_f}}. \quad (20) $$

The filter should be designed by respecting the following [3]:

$$ N_f > [\max(\beta_M, \alpha_N)], \quad (21) $$

where $[\cdot]$ is the ceiling operator. It has been shown in [5] that the cut-off frequency $\omega_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 [33]. This means that a priori estimation or tuning of the system bandwidth may be required to get better results.

The error function $\epsilon(t)$ is now replaced by its filtered version $\epsilon_f(t)$:

$$ \epsilon_f(t) = y_f^*(t) - \phi_f^*(t)^T \hat{\theta} \quad (22) $$

where the filtered regression vector $\phi_f^*(t)$ is

$$ \phi_f^*(t)^T = \begin{bmatrix} p^{\beta_0} u_f(t), \dots, p^{\beta_M} u_f(t), \\ -p^{\alpha_1} y_f^*(t), \dots, -p^{\alpha_N} y_f^*(t) \end{bmatrix}. \quad (23) $$

and the filtered input $u_f(t)$ and output $y_f(t)$ are defined as:

$$ \begin{aligned} u_f(t) &= F(p)u(t) \\ y_f^*(t) &= F(p)y^*(t). \end{aligned} \quad (24) $$

Inspired by the recursive least-squares [19], the RLSSVF algorithm is obtained for fractional order systems:

$$ \left\{ \begin{aligned} \hat{\theta}(k) &= \hat{\theta}(k-1) + L(k) \left[ y_f^*(kT_s) - \phi_f^*(kT_s)^T \hat{\theta}(k-1) \right] \\ L(k) &= \frac{F(k-1)\phi_f^*(kT_s)}{\lambda(k) + \phi_f^*(kT_s)^T F(k-1)\phi_f^*(kT_s)} \\ F(k) &= \frac{1}{\gamma(k)} \left[ I - L(k)\phi_f^*(kT_s)^T \right] F(k-1) \end{aligned} \right. \quad (25) $$

where $\gamma$ 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 $\gamma$ may actually be a function of time and be adaptive as identification is running [40]. In this study, it will be taken as a constant and equal to unity.

#### 3.2 Recursive Prediction Error Method (RPEM)

By following the prediction error method, a suitable error function is given by the output error:

$$ \epsilon(t) = y(t) - \hat{y}(t), \quad (26) $$

where the estimated output $\hat{y}(t)$ is computed as:

$$ \hat{y}(t) = G(p, \hat{\theta})u(t). \quad (27) $$

The quadratic error function defined in (18) is minimized. By taking the gradient of the criterion:

$$ \frac{\partial J(\hat{\theta})}{\partial \hat{\theta}} = \sum_{k=1}^{t} \psi_{\rho}(kT_s) \epsilon(kT_s) \quad (28) $$
with

$$ \psi_{\rho}(kT_s) = \frac{\partial \epsilon(kT_s)}{\partial \theta}, \qquad (29) $$

its error sensitivity is computed by:

$$ \frac{\partial \epsilon(kT_s, \hat{\theta})}{\partial \hat{\theta}} = - \left[ \frac{\partial \hat{y}}{\partial b_0}, \dots, \frac{\partial \hat{y}}{\partial b_M}, \frac{\partial \hat{y}}{\partial a_1}, \dots, \frac{\partial \hat{y}}{\partial a_N} \right]^T \bigg|_{t=kT_s} \quad (30) $$

where

$$ \frac{\partial \hat{y}}{\partial b_i} = \frac{p^{\beta_i}}{1 + \sum_{j=1}^{N} a_j p^{\alpha_j}} u(kT_s) \quad (31) $$

and

$$ \frac{\partial \hat{y}}{\partial a_{\ell}} = - \frac{\sum_{i=0}^{M} b_i p^{\beta_i + \alpha_{\ell}}}{\left(1 + \sum_{j=1}^{N} a_j p^{\alpha_j}\right)^2} u(kT_s). \quad (32) $$

for $i = 0, \dots, M$ and $\ell = 1, \dots, N$.

Inspired by the recursive least-squares [19], the RPEM algorithm is obtained for fractional order systems:

$$ \left\{ \begin{array}{l} \epsilon(kT_s) = y^*(kT_s) - \hat{y}(kT_s) \\ \hat{\theta}(k) = \hat{\theta}(k-1) + \gamma_{\rho} R^{-1}(k) \psi_{\rho}(kT_s) \epsilon(kT_s) \\ R(k) = R(k-1) + \gamma_{\rho} [\psi_{\rho}(kT_s) \psi_{\rho}^{T}(kT_s) - R(k-1)] \end{array} \right. \quad (33) $$

where $\gamma_{\rho}$ is a refining gain analogous to the step in a gradient descent. This forgetting factor $\gamma_{\rho}$ may actually be adapted during the identification or be a function of time. However, the determination of an adequate $\gamma_{\rho}$ may prove to be non-trivial, as stated in [18]. For the present work, it will be taken as a constant.

### 3.3 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 derivative is needed, therefore an extended error with long memory is more suited:

$$ \tilde{\epsilon}(kT_s) = [\epsilon(0), \epsilon(T_s), \epsilon(2T_s), \dots, \epsilon(kT_s)]^T. \quad (34) $$

This error signal $\tilde{\epsilon}$ 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 $\tilde{\psi}_{\rho}(k, \hat{\theta})$ is now a matrix:

$$ \tilde{\psi}_{\rho}(kT_s) = \begin{bmatrix} \frac{\partial \epsilon(0)}{\partial b_0} & \frac{\partial \epsilon(T_s)}{\partial b_0} & \frac{\partial \epsilon(2T_s)}{\partial b_0} & \dots & \frac{\partial \epsilon(kT_s)}{\partial b_0} \\ \frac{\partial \epsilon(0)}{\partial b_1} & \frac{\partial \epsilon(T_s)}{\partial b_1} & \frac{\partial \epsilon(2T_s)}{\partial b_1} & \dots & \frac{\partial \epsilon(kT_s)}{\partial b_1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \epsilon(0)}{\partial a_N} & \frac{\partial \epsilon(T_s)}{\partial a_N} & \frac{\partial \epsilon(2T_s)}{\partial a_N} & \dots & \frac{\partial \epsilon(kT_s)}{\partial a_N} \end{bmatrix} \quad (35) $$

where $k$ indicates the present iteration.

By introducing the extended measured output $\tilde{Y}^*(k)$ and the extended estimated output $\tilde{Y}(k)$ as:

$$ \tilde{Y}^*(kT_s) = [y^*(0), y^*(T_s), y^*(2T_s), \dots, y^*(kT_s)]^T \quad (36) $$

$$ \tilde{Y}(kT_s) = [\hat{y}(0), \hat{y}(T_s), \hat{y}(2T_s), \dots, \hat{y}(kT_s)]^T, \quad (37) $$

the long memory recursive prediction-error method (LM-RPEM) is proposed for fractional order systems:

$$ \left\{ \begin{array}{l} \tilde{\epsilon}(kT_s) = \tilde{Y}^*(kT_s) - \tilde{Y}^*(kT_s) \\ \hat{\theta}(k) = \hat{\theta}(k-1) + \gamma_{\rho} R^{-1}(k) \tilde{\psi}_{\rho}(kT_s) \tilde{\epsilon}(kT_s) \\ R(k) = R(k-1) + \gamma_{\rho} [\tilde{\psi}_{\rho}(kT_s) \tilde{\psi}_{\rho}^{T}(kT_s) - R(k-1)] \end{array} \right. \quad (38) $$

where again, $\gamma_{\rho}$ is a refining gain analogous to the step in a gradient descent. As previously stated in section 3.2, $\gamma_{\rho}$ will be taken as a constant.

### 3.4 Simulation results for recursive coefficient estimation

#### 3.4.1 True system description

Let us consider the following true system model:

$$ G(p) = \frac{K}{1 + 2\zeta \left(\frac{p}{\omega_0}\right)^{\nu} + \left(\frac{p}{\omega_0}\right)^{2\nu}} = \frac{1}{1 - 1.414p^{0.5} + 2p^1} \quad (39) $$

with $K = 1$, $\zeta = -0.5$, $\omega_0 = 0.5\text{rad/s}$ and $\nu = 0.5$. Such a general structure represents a second kind fractional order system [23], and has a resonant peak around $\omega_r = 0.45\text{rad/s}$ (see Fig. 1).

The model structure is chosen such as the true one (39) and is expressed by:

$$ G(p) = \frac{b_0}{1 + a_1 p^{\nu} + a_2 p^{2\nu}}, \quad (40) $$

with differentiation orders assumed known and fixed just like the true system (39) ($\nu = 0.5$). The aim is to estimate the coefficients of the system and to see the statistical properties of the proposed LMRPEM algorithm
<figure><img src="image_3.png" /><figcaption>Fig. 1 Bode diagram of the true system (39)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2 Input/output data used for recursive system identification</figcaption></figure>

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.01$s. All simulations are initialized with $\theta^T = [1, 0.5, 2.5]$ and $\nu = 0.5$ such as the true system (39). Fig. 2 shows the input used for simulations, as well as the noise-free output and the noisy output with a SNR = 15dB.

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<fnref n="2" />.

<table><caption>Table 1 Monte Carlo simulation results (100 runs) for coefficient estimation with different SNR. $\tilde{\theta}$ denotes the parameter mean value and $\hat{\sigma}_{\theta}$ the standard deviations</caption><thead><tr><th rowspan="2">True<br>Method</th><th rowspan="2">SNR</th><th colspan="2">b<sub>0</sub></th><th colspan="2">a<sub>2</sub></th><th colspan="2">a<sub>1</sub></th></tr><tr><th>̂b<sub>0</sub></th><th>̂σ<sub>b<sub>0</sub></sub></th><th>̂a<sub>2</sub></th><th>̂σ<sub>a<sub>2</sub></sub></th><th>̂a<sub>1</sub></th><th>̂σ<sub>a<sub>1</sub></sub></th></tr></thead><tbody><tr><td rowspan="3">RLSSVF</td><td>20</td><td>0.866</td><td>0.014</td><td>1.740</td><td>0.017</td><td>-1.317</td><td>0.016</td></tr><tr><td>15</td><td>0.644</td><td>0.026</td><td>1.357</td><td>0.030</td><td>-1.110</td><td>0.059</td></tr><tr><td>10</td><td>0.295</td><td>0.035</td><td>0.819</td><td>0.042</td><td>-0.764</td><td>0.138</td></tr><tr><td rowspan="3">RPEM</td><td>20</td><td>1.001</td><td>0.010</td><td>1.999</td><td>0.034</td><td>-1.414</td><td>0.031</td></tr><tr><td>15</td><td>0.999</td><td>0.020</td><td>1.985</td><td>0.051</td><td>-1.403</td><td>0.050</td></tr><tr><td>10</td><td>1.000</td><td>0.035</td><td>1.954</td><td>0.098</td><td>-1.374</td><td>0.102</td></tr><tr><td rowspan="3">LMRPEM</td><td>20</td><td>0.999</td><td>0.004</td><td>1.999</td><td>0.006</td><td>-1.415</td><td>0.006</td></tr><tr><td>15</td><td>0.999</td><td>0.009</td><td>2.001</td><td>0.011</td><td>-1.416</td><td>0.012</td></tr><tr><td>10</td><td>0.998</td><td>0.016</td><td>2.001</td><td>0.019</td><td>-1.417</td><td>0.022</td></tr></tbody></table>

### 3.4.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 cut-off frequency of $\omega_f = 10$rad/s and as $\max(\beta_M, \alpha_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 $\omega_f$ and the order $N_f$ will be required.

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, the greater the estimation bias and the greater the standard deviations. 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 cut-off frequency of the filter, however, it should be noted that a frequency $\omega_f$ extremely close to the system bandwidth would deteriorate the estimation in a practical scenario.

The RPEM method ($\gamma_\rho$ 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 = 20dB and 10dB.

The LMRPEM method ($\gamma_\rho$ was set to 0.01) has the most consistent estimation as the estimation has no bias and the estimation variances are the lowest whatever the noise level. This method has significantly reduced the parameter deviations.

<footnote n="2">² The new version of the CRONE toolbox is an object-oriented version with several classes defined for fractional models (LTI, explicit form, implicit form, ZPK, state-space representation, etc.). This CRONE toolbox is freely available on http://archive.ims-bordeaux.fr/CRONE/toolbox/</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 3 Criterion J for different commensurate orders $\nu$ ($SNR = 20$ dB)</figcaption></figure>

## 3.5 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 $l_2$-norm of the normalized output error:

$$J(\hat{\theta}) = 10\log \frac{\|\tilde{\epsilon}(NT_s, \hat{\theta})\|^2}{\|\tilde{Y}(NT_s), \hat{\theta}\|^2} \quad (41)$$

where the long memory error $\tilde{\epsilon}(NT_s, \hat{\theta})$ is defined as in (34):

$$\tilde{\epsilon}(NT_s, \hat{\theta}) = \tilde{Y}^*(NT_s) - \tilde{Y}(NT_s, \hat{\theta}) \quad (42)$$

are given by 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 NSR = -20 dB. As expected, the criterion minimum is obtained at $\nu = 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 $\nu = 0.5$ and equals -20dB which exactly corresponds to the NSR. Other values of $\nu$ lead to non-negligible modeling errors as the coefficient estimation would not converge to the true parameters. For example, at $\nu = 1$, the cost function is around -10.01dB, which corresponds to a modeling error around 10dB.

This simulation result motivates to estimate the fractional differentiation orders, as they may considerably influence the estimation.

## 4 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 previous section for coefficient estimation, it is proposed to use this method for the differentiation order estimation. Then, two variants are proposed to identify all parameters, both the coefficients and the differentiation orders in a two stage algorithm.

### 4.1 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:

$$\epsilon(t) = y^*(t) - \hat{y}(t), \quad (43)$$

and differentiation orders may be recursively adjusted by adding a correction term to the previous iteration:

$$\mu(k+1) = \mu(k) + \gamma_{\mu} R_{\mu}^{-1} \psi_{\mu}(kT_s) \epsilon(kT_s) \quad (44)$$

where:

$$\psi_{\mu}(kT_s) = \frac{\partial \epsilon(kT_s)}{\partial \mu} = -\frac{\partial \hat{y}}{\partial \mu}. \quad (45)$$

Recall that the LMRPEM method has given the best results for estimating the coefficients for fractional order systems. Consequently, an extended error is used to take into account the long memory behavior of the system:

$$\tilde{\epsilon}(kT_s) = [\epsilon(0) \epsilon(T_s) \epsilon(2T_s) \dots \epsilon(kT_s)]^T, \quad (46)$$

which turns the error sensitivity into a matrix:

$$\tilde{\psi}_{\mu}(kT_s) = \begin{bmatrix} \frac{\partial \epsilon(0)}{\partial \beta_0} & \frac{\partial \epsilon(T_s)}{\partial \beta_0} & \frac{\partial \epsilon(2T_s)}{\partial \beta_0} & \dots & \frac{\partial \epsilon(kT_s)}{\partial \beta_0} \\ \frac{\partial \epsilon(0)}{\partial \beta_1} & \frac{\partial \epsilon(T_s)}{\partial \beta_1} & \frac{\partial \epsilon(2T_s)}{\partial \beta_1} & \dots & \frac{\partial \epsilon(kT_s)}{\partial \beta_1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \epsilon(0)}{\partial \alpha_N} & \frac{\partial \epsilon(T_s)}{\partial \alpha_N} & \frac{\partial \epsilon(2T_s)}{\partial \alpha_N} & \dots & \frac{\partial \epsilon(kT_s)}{\partial \alpha_N} \end{bmatrix}. \quad (47)$$

For commensurate order estimation, the gradient boils down to:

$$\tilde{\psi}_{\mu}(kT_s) = \tilde{\psi}_{\nu}(kT_s) = \begin{bmatrix} \frac{\partial \epsilon(0)}{\partial \nu} & \frac{\partial \epsilon(T_s)}{\partial \nu} & \frac{\partial \epsilon(2T_s)}{\partial \nu} & \dots & \frac{\partial \epsilon(kT_s)}{\partial \nu} \end{bmatrix}. \quad (48)$$
The commensurate order output sensitivity $\frac{\partial \hat{y}}{\partial \nu}$ is expressed by:

$$ \frac{\partial \hat{y}}{\partial \nu} = \ln(p) \frac{\sum_{i=0}^{M} ib_i p^{i\nu} \sum_{i=0}^{M} \sum_{j=1}^{N} (i-j)b_i a_j p^{(i+j)\nu}}{\left(1 + \sum_{j=1}^{N} a_i p^{\alpha_i}\right)^2} u(t). \quad (49) $$

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:

$$ \left\{ \begin{array}{l} \tilde{\epsilon}(kT_s) = \tilde{Y}^*(kT_s) - \tilde{Y}^*(kT_s) \\ \hat{\nu}(k) = \hat{\nu}(k-1) + \gamma_{\nu} R^{-1}(k) \tilde{\psi}_{\nu}(kT_s) \tilde{\epsilon}(kT_s) \\ R(k) = R(k-1) + \gamma_{\nu} \left[ \tilde{\psi}_{\nu}(kT_s) \tilde{\psi}_{\nu}^{T}(kT_s) - R(k-1) \right] \end{array} \right. \quad (50) $$

## 4.2 All parameter system identification: hybrid algorithms

The best previous method for coefficient estimation is combined with differentiation order estimation for all parameter estimation of a fractional order model. This hybrid LMRPEM algorithm has a coefficient estimation step and a differentiation order estimation step.

<caption>**Algorithm 1: Hybrid LMRPEM algorithm**</caption>

Step 1: Initialization
$k = 0$
Initialize parameter vector $\theta(0) = [\rho(0)^T \nu(0)^T]^T$
Step 2: Recursive algorithm
for $k = 1$ : end

(a) Coefficient estimation
Assuming the order $\nu(k-1)$ known,
compute $\rho(k)$ from equation (38)

(b) Differentiation order estimation
Assuming the coefficient $\rho(k)$ known,
compute $\nu(k)$ from equation (50)

As stated by Young in [44], the tuning of the forgetting factor $\gamma$ in recursive PEM methods is non-trivial and influences the convergence of the parameters. Consequently, algorithm 1 requires the tuning of two forgetting factors (one for the coefficients and one for the differentiation orders).

## 4.3 Simulation results for recursive all parameter estimation

The same system (39) as detailed in section 3.4.1 is used with different noise levels: $\infty$ without noise, 20dB,

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Parameter estimation evolution with hybrid LMR-PEM alg. 1 for SNR = 20dB</figcaption></figure>

and 15dB. The tuning of the $\gamma$ gains was performed empirically and is given in Table 2.

<caption>**Table 2** Hybrid LMRPEM algorithm parameters</caption>

Hybrid algorithm 1 $\gamma_{\rho} = 0.005$ and $\gamma_{\nu} = 0.05$

The parameter evolution for an SNR = 20dB are plot in Fig. 4 for Alg. 1. On the commensurate order estimation, hybrid LMRPEM Alg. 1 exhibits some oscillations at the beginning. However, the hybrid LMR-PEM algorithm tends to the true commensurate order $\nu = 0.5$.

Monte Carlo Simulation analysis are 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 Alg. 1 has converged to the true commensurate order, and therefore the coefficient estimations is pretty accurate in mean. As the noise level increases, Hybrid LMRPEM Alg. 1 presents consistent estimates, with minimum 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 derivative 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 computation time by the end of the simulation, which would
<table><caption>**Table 3** Monte Carlo simulation results (100 runs) for all parameter estimation with different SNR. $\tilde{\theta}$ denotes the parameter mean value and $\hat{\theta}_{\theta}$ the standard deviations</caption><thead><tr><th></th><th>b<sub>0</sub></th><th>a<sub>2</sub></th><th>a<sub>1</sub></th><th>ν</th></tr></thead><tbody><tr><th>True</th><td>1.000</td><td>2.000</td><td>-1.414</td><td>0.5</td></tr><tr><th>SNR</th><td>0.991</td><td>2.001</td><td>-1.431</td><td>0.497</td></tr><tr><th>Hybrid LMRPEM</th><td>∞</td><td>0.987</td><td>0.023</td><td>0.496</td></tr><tr><th>Alg. 1</th><td>20</td><td>0.989</td><td>0.029</td><td>0.007</td></tr><tr><th></th><td>15</td><td>0.989</td><td>0.029</td><td>0.009</td></tr></tbody></table>

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 [35].

## 5 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 was previously proposed in [42] 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 figure 5). All model developments are detailed in [7].

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Thermal two-port network</figcaption></figure>

The bronchus is modeled as an air cylinder with an intermediate length of $L = 0.0236\text{m}$ [27] and radius $r \approx 1\text{mm}$ [26]. Impedance approximations are given by:

$$Z_{1}(s) = Z_{2}(s) \approx \frac{150240}{1 + 2.5158\sqrt{s}} \quad (51)$$

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Input and output data for bronchus application (zoom on 500s of the 2000s dataset)</figcaption></figure>

and

$$Z_{3}(s) \approx \frac{11869}{s} \quad (52)$$

For an insulated output case ($\phi_{out} = 0$), the transfer function relating input and output temperatures is given by:

$$G_{\text{bronchus}}(s) = \frac{2.516\sqrt{s} + 1}{12.66s + 2.516\sqrt{s} + 1} \quad (53)$$

The input signal is a pseudo random binary sequence oscillating between -1 and 1 and containing 2000 values with a sampling time $T_s = 1\text{s}$. The initialization is $\theta^T = [2, 5, 5, 10]$ and $\nu = 0.35$. The noise was set to SNR = 20 dB. The input/output data is shown in figure 6.

This model is set as the true system (53):

$$G(s) = \frac{b_1 s^{\nu} + b_0}{a_2 s^{2\nu} + a_1 s^{\nu} + 1} \quad (54)$$

Hybrid LMRPEM Alg. 1 is used to estimate the parameters as it is the most efficient algorithm for all parameter estimation (see section 4.3). $\gamma_{\rho}$ and $\gamma_{\nu}$ are the same as given in Table 2. The recursive coefficient and differentiation order estimations are shown in figure 7. The commensurate order convergence is slow but very accurate towards the system true commensurate order $\nu = 0.5$. The coefficient convergence exhibit fairly
<figure><img src="image_9.png" /><figcaption>Fig. 7 Parameter estimation for bronchus application</figcaption></figure>

slow convergence rate, especially for $a_2$, however as it tends to its true value towards the end of the simulation. Therefore, the longer the dataset, the better the convergence to the true parameters.

Finally, at the last iteration, the estimated model is expressed as the following:

$$ \hat{G}(s) = \frac{2.519s^{0.502} + 1.013}{12.66s^{1.005} + 2.601s^{0.502} + 1}, \quad (55) $$

an estimated model that has well converged to the true parameters such as predicted by hybrid LMRPEM Alg. 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.

## 6 Conclusions and perspectives

A new method called LMRPEM is proposed for parameter recursive system identification of continuous-time fractional order models. Least squares (RLSSVF), 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 method have proven to produce the lesser parameter variance as compared to other methods.

Then, for all parameter estimation, especially for differentiation order estimation, a new hybrid algorithm is introduced: the hybrid LMRPEM for both coefficient and differentiation order estimation. The LMRPEM 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 behaviour.

Finally, a biological example was proposed: heat transfer in lungs is modeled on a human main bronchus by using thermal two-port network formalism with fractional models. Applying the hybrid LMRPEM method have provided consistent estimates without bias and very low variance.

Research perspectives may include further studies regarding a method to tune the forgetting factor $\gamma$ 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 sheeps will be carried out at the Bordeaux University Hospital and IHU Liryc Institute for online system identification of lung heat transfers.

## 7 Declarations

Funding: Not applicable

Conflicts of interest/Competing interests: The authors have no relevant financial or non-financial interests to disclose.

Code availability (software application or custom code): Not applicable

Authors' Contributions - in relation with your submission: J.-F. Duhé has developed the recursive algorithms with S. Victor. He also made the simulations. S. Victor, P. Melchior, Y. Abdelmoumen and F. Roubertie have contributed to the biological and physiological results.

Ethics approval (include appropriate approvals or waivers): Not applicable

Consent to participate (include appropriate statements): Not applicable

Consent for publication (include appropriate statements): Not applicable

The manuscript has no associated data.

## References

1. Battaglia, J.L., Le Lay, L., Batsale, J.C., Oustaloup, A., Cois, O.: Heat flux estimation through inverted non integer identification models. International Journal of Thermal Science 39(3), 374-389 (2000). DOI 10.1016/S1290-0729(00)00220-9
2. Chen, F., Garnier, H., Gilson, M., Zhuan, X.: Frequency domain identification of continuous-time output-error models with time-delay from relay feedback tests. Automatica **98**, 180–189 (2018)

3. Cois, O., Oustaloup, A., Poinot, T., Battaglia, J.L.: Fractional state variable filter for system identification by fractional model. In: 6th European Control Conference ECC’01. Porto, Portugal (2001)

4. Das, S., Sivaramakrishna, M., Das, S., Biswas K. and Goswami, B.: Characterization of a fractional order element realized by dipping a capacitive type probe in polarizable medium. In: Symposium on Fractional Signals and Systems. Lisbon, Portugal (2009)

5. De Wit, C.: Recursive estimation of the continuous-time process parameters. In: 1986 25th IEEE Conference on Decision and Control, pp. 2016–2020 (1986). DOI 10.1109/CDC.1986.267390

6. Djouambi, A., Besançon, A.V., Charef, A.: Fractional system identification using recursive algorithms approach. In: 2007 European Control Conference (ECC), pp. 1436–1441 (2007). DOI 10.23919/ECC.2007.7068707

7. Duhe, J., Victor, S., Melchior, P., Abdelmoumen, Y., Roubertie, F.: Modeling thermal systems with fractional models: human bronchus application. Nonlinear Dynamics accepted (2022)

8. Eddine, A., Huard, B., Gabano, J.D., Poinot, T.: Initialization of a fractional order identification algorithm applied for lithium-ion battery modeling in time domain. Commun Nonlinear Sci Numer Simulat **59**, 375–386 (2018)

9. Elwakil, A.: Fractional-order circuits and systems: An emerging interdisciplinary research area. IEEE Circuits and Systems Magazine **10**(4), 40–50 (2010). DOI 10.1109/MCAS.2010.938637

10. Garnier, H., Wang, L.: Identification of continuous-time models from sampled data. Springer-Verlag (2008)

11. Garrappa, R., Kaslik, E., Popolizio, M.: Evaluation of fractional integrals and derivatives of elementary functions: Overview and tutorial. Mathematics **7**(407) (2019). DOI 10.3390/math7050407

12. Grünwald, A.: Uber begrenzte Derivationen und deren Anwendung. Zeitschrift für Mathematik und Physik pp. 441–480 (1867)

13. Idiou, D., Charef, A., Djouambi, A., Voda, A.: Parameters and order identification of the fundamental linear fractional systems of commensurate order. In: The Second International Conference on Electrical Engineering and Control Application, ICEECA 2014. Constantine, Algeria (2014)

14. Ionescu, C., Copot, D., De Keyser, R.: Respiratory impedance model with lumped fractional order diffusion compartment. IFAC Proceedings Volumes **46**(1), 260 – 265 (2013). DOI 10.3182/20130204-3-FR-4032.00084. 6th IFAC Workshop on Fractional Differentiation and Its Applications

15. Kilbas, A.A., Srivastava, H.M., Trujillo, J.J.: Theory and Applications of Fractional Differential Equations, Volume 204 (North-Holland Mathematics Studies). Elsevier Science Inc., New York, NY, USA (2006)

16. Krishna, B.: Studies on fractional order differentiators and integrators: A survey. Signal Processing **91**(3), 386 – 426 (2011). DOI 10.1016/j sigpro.2010.06.022

17. Letnikov, A.: Theory of differentiation of arbitrary order (russian). Matematiceskij Sbornik (Moscou) **3**(1), 1–68 (1868)

18. Ljung, L.: Analysis of a general recursive prediction error identification algorithm. Automatica **17**(1), 89–99 (1981). DOI 10.1016/0005-1098(81)90086-8

19. Ljung, L.: System identification – Theory for the user, 2 edn. Prentice-Hall, Upper Saddle River, N.J., USA (1999)

20. Ljung, L., Chen, T., Mu, B.: A shift in paradigm for system identification. International Journal of Control **93**(2), 173–180 (2020). DOI 10.1080/00207179.2019.1578407

21. Magin, R., Ovadia, M.: Modeling the cardiac tissue electrode interface using fractional calculus. 2nd IFAC Workshop on Fractional Differentiation and its Applications **39**(11), 302 – 307 (2006). DOI 10.3182/20060719-3-PT-4902.00056

22. Maillet, D., André, S., Batsale, J., Degiovanni, A., Moyne, C.: Thermal Quadrupoles: Solving the Heat Equation through Integral Transforms. Loyola Symposium Series. John Wiley & Sons (2000)

23. Malti, R., Moreau, X., Khemane, F., Oustaloup, A.: Stability and resonance conditions of elementary fractional transfer functions. Automatica **47**(11), 2462–2467 (2011). DOI 10.1016/j.automatica.2011.08.029

24. Malti, R., Sabatier, J., Akçay, H.: Thermal modeling and identification of an aluminium rod using fractional calculus. In: 15th IFAC Symposium on System Identification (SYSID’2009), pp. 958–963. St Malo, France (2009). DOI 10.3182/20090706-3-FR-2004.00159

25. Matignon, D.: Stability properties for generalized fractional differential systems. ESAIM proceedings - Systèmes Différentiels Fractionnaires - Modèles, Méthodes et Applications **5** (1998)

26. McFawn, P., Mitchell, H.: Bronchial compliance and wall structure during development of the immature human and pig lung. European Respiratory Journal **10**(1), 27–34 (1997)

27. Mi, W., Zhang, C., Wang, H., Cao, J., Li, C., Yang, L., Guo, F., Wang, X., Yang, T.: Measurement and analysis of the tracheobronchial tree in chinese population using computed tomography. PLOS ONE **10**(4), 1–14 (2015). DOI 10.1371/journal.pone.0123177

28. Moze, M., Sabatier, J.: LMI tools for stability analysis of fractional systems. In: 20th ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, IDETC/CIE’05, pp. 1–9. Long Beach, CA (2005)

29. Nakagawa, M., Sorimachi, K.: Basic characteristics of a fractance device. IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences **75**, 1814–1819 (1992)

30. Narang, A., Shah, S., Chen, T.: Continuous-time model identification of fractional-order models with time delays. IET Control Theory & Applications **5**(7), 900–912 (2011). DOI 10.1049/iet-cta.2010.0718

31. Oldham, K., Spanier, J.: The replacement of Fick’s laws by a formulation involving semidifferentiation. Journal of Electroanalytical Chemistry and Interfacial Electrochemistry **26**(2-3), 331–341 (1970). DOI 10.1016/S0022-0728(70)80316-3

32. Oldham, K., Spanier, J.: The fractional calculus - Theory and Applications of Differentiation and Integration to Arbitrary Order. Academic Press, New-York and London (1974)

33. Padilla, A.: Recursive identification of continuous-time systems with time-varying parameters. Theses, Université de Lorraine (2017)

34. Pillonetto, G., Dinuzzo, F., Chen, T., De Nicolao, G., Ljung, L.: Kernel methods in system identification, machine learning and function estimation: A survey. Automatica **50**(3), 657–682 (2014). DOI https://doi.org/10.1016/j.automatica.2014.01.001
35. Podlubny, I.: Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Academic Press, San Diego (1999)

36. Rivero, M., Rogosin, S., Tenreiro Machado, J., Trujillo, J.: Stability of fractional order systems. Mathematical Problems in Engineering **2013**, 356215 (2013). DOI 10.1155/2013/356215

37. Rodrigues, S., Munichandraiah, N., Shukla, A.K.: A review of state of charge indication of batteries by means of A.C. impedance measurements. Journal of Power Sources **87**(1-2), 12–20 (2000). DOI 10.1016/S0378-7753(99)00351-1

38. Samko, S., Kilbas, A., Marichev, O.: Fractional integrals and derivatives: theory and applications. Gordon and Breach Science (1993)

39. Sommacal, L., Melchior, P., Oustaloup, A., Cabelguen, J.M., Ijspeert, A.: Fractional multi-models of the frog gastrocnemius muscle. Journal of Vibration and Control **14**(9-10), 1415–1430 (2008). DOI 10.1177/1077546307087440

40. Sun, X., Ji, J., Ren, B., Xie, C., Yan, D.: Adaptive forgetting factor recursive least square algorithm for online identification of equivalent circuit model parameters of a lithium-ion battery. Energies **12**(12) (2019). DOI 10.3390/en12122242

41. Victor, S., Malti, R., Garnier, H., Oustaloup, A.: Parameter and differentiation order estimation in fractional models. Automatica **49**(4), 926–935 (2013). DOI 10.1016/j.automatica.2013.01.026

42. Victor, S., Melchior, P., Pellet, M., Oustaloup, A.: Lung thermal transfer system identification with fractional models. IEEE Transactions on Control Systems Technology **28**(1), 172–182 (2020). DOI 10.1109/TCST.2018.2877606

43. Wang, L., Zhao, W.: System identification: New paradigms, challenges, and opportunities. Acta Automatica Sinica **39**(7), 933–942 (2013). DOI https://doi.org/10.1016/S1874-1029(13)60062-2

44. Young, P.: Parameter estimation for continuous-time models – a survey. Automatica **17**(1), 23–29 (1981)

45. Zeng, C., Liang, S.: Comparative study of discretization zero dynamics behaviors in two multirate cases. International Journal of Control, Automation and Systems **13**(4), 831–842 (2015). DOI 10.1007/s12555-014-0115-3