# Recursive System Identification of Continuous-Time Fractional Systems for All Parameter Estimation

Jean-François Duhé, Stéphane Victor, Pierre Melchior, Youssef

Abdelmoumen, François Roubertie

► To cite this version:

Jean-François Duhé, Stéphane Victor, Pierre Melchior, Youssef Abdelmoumen, François Rouber-tie. Recursive System Identification of Continuous-Time Fractional Systems for All Parameter Estimation. IEEE Transactions on Control Systems Technology, 2024, 32 (6), pp.2037 - 2049. 10.1109/TCST.2024.3407580. [hal-04716247](https://hal.science/hal-04716247v1)

HAL Id: hal-04716247

https://hal.science/hal-04716247v1

Submitted on 12 Jan 2026

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Recursive system identification of
continuous-time fractional systems for all
parameter estimation

Jean-François Duhé, Stéphane Victor, Pierre Melchior, Youssef Abdelmoumen and François Roubertie

Abstract

Fractional order systems have proven to be useful to well model diffusion or propagation phenomena. Recursive or online system identification of continuous-time fractional models is explored in this paper. When differentiation orders are known, only the coefficients are to be estimated: the classic recursive methods of least-squares, prediction error method and instrumental variable are adapted for fractional models. They are then compared to our new long memory prediction error method to prove its efficiency. When differentiation orders are unknown, which is often the case in practice, two-stage algorithms are proposed for both coefficient and differentiation order estimation. A single method of differentiation order estimation is proposed, which is then combined with the two best coefficient estimation methods (Long Memory Prediction Error Method and Instrumental Variable) to create two hybrid algorithms. The performances of these algorithms are compared through Monte-Carlo simulations in order to highlight the influence of the parameter estimation in a more complex scenario. Finally, recursive identification is applied to a simulation example of a thermal lung impedance.

Index Terms

Long Memory Prediction Error Method, Instrumental Variable, Recursive identification, Monte-Carlo simulation, continuous-time models, fractional calculus, fractional order model, thermal system, lung modeling

# I. INTRODUCTION

TRADITIONAL calculus definitions of integrals and derivatives can be generalized in order to conceive a general operator: the differ-integral (commonly called fractional order derivative). The creation of such an operator dates back to Riemann [1] and Liouville [2]. This concept remained limited to an abstract mathematical operator which contains the classic definitions of integrals and derivatives as well as a transition between such operators. However, during the 20th century this tool was explored as a way to model diffusive phenomena for complex geometries [3], [4]. Fractional order models have been shown to model an infinity of dynamics contained in a system, which can be observed by means of recursive cell association [5], [6]. This concept can be applied for dynamical analysis of fractal geometries, which have proven to be present in biology [7]. During the sixties, the study of mechanical and dielectric relaxation processes led to fractional order models in the complex plane [8]. Heat equation leads to half-order fractional integrator behavior for diffusion [9], [10]. Circuit models have long been used as a tool for creating lumped-parameter models [11]. Fractional integrators may be seen as an element between impedance and capacitance called constant-phase element [12]–[14]. Constant-phase elements have shown to be useful in biological scenarios, such as intestine tissue dynamics [15], cardiac tissue [16] and lung mechanical impedance [17], [18]. They are also used in electro-chemistry, notably for battery impedance [19] and porous films. The enhanced flexibility in the frequency domain provided by these operators may also be exploited in filter-design such as the generalized fractional low-pass filter [20]. Artificial potential fields, both for attractive and repulsive potential fields, have also been explored for unmanned autonomous systems such as cars [21], [22] and drones [23].

The main goal of system identification is to give an accurate estimation of the model parameters used to represent a real system by analyzing experimental data. There is a vast literature concerning the analysis and identification of discrete-time models [24]. However, it should be noted that such type of models may lead to parameters that lack physical meaning. The sampling zero problem [25] related to discretization should neither be neglected. As a physically meaningful model estimation is desired, 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 gathered in the regression matrix $\phi$ [26]. If the models

"This work has no research funding"

J.-F. Duhé is with Universidad de Panamá, Panamá, Panamá (e-mail: jean-f.duhe-p@up.ac.pa).

S. Victor is with the Automatic Control group, CRONE team at the IMS laboratory, UMR 5218, CNRS and Université de Bordeaux, 33405 Talence cedex, France (e-mail: stephane.victor@ims-bordeaux.fr).

P. Melchior is with the Automatic Control group, CRONE team at the IMS laboratory, UMR 5218, CNRS and Université de Bordeaux, 33405 Talence cedex, France (e-mail: pierre.melchior@ims-bordeaux.fr)

Y. Abdelmoumen is with Clinique Saint-Augustin, 33000 Bordeaux Cedex, France (e-mail: dr.a.youssef@free.fr)

F. Roubertie is with IHU Liryc, Electrophysiology and Heart Modeling Institute, Fondation Bordeaux Université, 33000 Bordeaux Cedex, France and Centre de Recherche Cardio-Thoracique de Bordeaux, INSERM, 33000 Bordeaux Cedex, France (e-mail: rouberf@me.com)
to be identified are fractional order ones, this regression matrix will contain approximations of fractional order derivatives from experimental data.

Fractional order system identification began in the 1990s. The earliest proposed methods [27], [28] estimate coefficients of fractional order transfer functions. However, these first methods were limited to the Grünwald-Letnikov discrete-time definition of the fractional derivative. By taking continuous-time identification methods [29], methods relying on state-variable filters, least-squares optimization and instrumental variables have been introduced for single-input single-output (SISO) systems in [30], [31] and have been extended to multiple-inputs single-output systems in [32], [33]. These methods allow the use of different definitions of the operator and have been applied in different scenarios: identification of a thermal model on an aluminium rod [31], [34], [35], lithium-ion battery models [36] and even muscular relaxation in frogs [37], propagation of epidemics such as COVID-19 [38].

The presented methods are all proposed in off-line, which implies that the whole dataset of an experiment is used and optimized afterwards the data acquisition has been completed. However, it may be sought in some cases to perform an online estimation in which parameters will be estimated as each new measurement is acquired. One of the most famous online identification methods is the recursive least squares. Recursive prediction-error methods are also available in the literature (see e.g. [24]). By starting from Le Lay's method and online estimation concepts, a method to correctly estimate coefficients of a fractional order transfer function was proposed in [39]. [40] provides a recursive versions of LSSVF (least-squares with state variable filters) and RIVC (refined instrumental variable for continuous-time systems) for integer order (or rational) systems. Fractional order models may also require the estimation of one or several differentiation orders, which has already been studied in the off-line cases by means of gradient descent [31]. An extension of Djouambi's identification method has also been proposed to identify the fractional differentiation order by using a scheme similar to gradient descent [41].

The paper has two main objectives: coefficient and order estimation. In the first part, the main goal is to recursively estimate coefficients for fractional order continuous-time models. Well-known recursive identification methods (LSSVF, PEM (prediction error method) and optimal IV (instrumental variable)) are extended and adapted in order to deal with fractional order models. In addition to this, a new method derived from the prediction error method is proposed which is more suited to fractional order models. These four adapted methods will be tested by using Monte Carlo simulations on an academic simulation model in order to compare their performances. For the second part, the concept of prediction error method will be considered for the estimation of differentiation orders in order to propose an estimation scheme. Then, two hybrid algorithms are proposed to estimate both coefficients and differentiation orders. Coefficient estimation stage will be different for both proposed algorithms whereas the differentiation order estimation will be carried out by our proposed LMRPEM method (long memory prediction error method). They will also be tested by means of Monte Carlo simulations.

Besides the mathematical contributions, a thermal application is also proposed. Equivalent thermal impedance of a lung is considered as the plant to be identified and the methods proposed in this work are applied to estimate both coefficients and differentiation order. Indeed such a robust online system identification algorithms are necessary in the case of online identification of heat transfers in lungs during cardiac surgery such as proposed in [42] for offline system identification of black-box systems.

The paper is organized as follows. Section II introduces the basics of fractional calculus and fractional order transfer functions, Section III presents recursive fractional system identification methods for coefficient estimation. Then, Section IV introduces an estimation method for unknown fractional differentiation orders as long as coefficients for online parameter identification. Section V presents simulation results for an academic scenario. The lung thermal impedance application is presented in Section VI. Finally, conclusions are drawn in Section VII.

## II. FRACTIONAL CALCULUS CONTEXT

### A. Mathematical background

a) Fractional derivative definition: Throughout history, there have been different definitions for the fractional derivative [43]–[45]. The most well-known definition in practice is the Grünwald-Letnikov definition [46] and [47].

Let $a \in \mathbb{R}$ and a function $f \in C^{\infty}([a, +\infty))$, the set of infinitely continuously differentiable functions from $[a, +\infty)$ to $\mathbb{R}$. The fractional differentiation to an arbitrary order $\alpha \in \mathbb{R}^+$ of $f(t)$, in the sense of Grünwald-Letnikov, is defined by:

$$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 $p$ denotes the differential operator ($p = \frac{d}{dt}$), $\lfloor \cdot \rfloor$ stands for the floor operator and $\binom{\alpha}{j}$ stands for Newton's binomial coefficient generalized to real numbers:

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

where $\Gamma$ is Euler's gamma function defined as:

$$\Gamma(x) = \int_{0}^{\infty} e^{-t} t^{x-1} dt \quad (3)$$
for $x \in \mathbb{R} - \mathbb{Z}^{-}$.

For simplicity sake, let us assume that $a = 0$, which is the time from which the signal is acquired.

The fractional derivative (1) is usually the most used one as, first of all, the initial conditions of the operator are well- considered and secondly, it can easily be approximated by a discrete definition by replacing $h$ parameter for a sampling time $T_s$ and dropping the limit:

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

As stated in [48, section 7.4], the error terms are proportional to the sampling period. Therefore, to make the approximation error negligible, the sampling time must be sufficiently small<fnref n="1" />.

Note that system identification algorithms proposed in this paper are meant to be independent from time-domain simulation algorithms. Nevertheless, fractional order models should be correctly simulated with negligible simulation errors so that the parameter estimation are consistent.

For null initial conditions, the Laplace transform of the fractional derivative leads to a simple expression [48]:

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

b) Fractional order dynamic models: A fractional SISO model may relate its output $y(t)$ to input $u(t)$ through a fractional order differential equation:

$$\begin{aligned} y(t) + a_1 p^{\alpha_1} y(t) + \dots + a_{m_A} p^{\alpha_{m_A}} y(t) = & \\ b_0 p^{\beta_0} u(t) + b_1 p^{\beta_1} u(t) + \dots + b_{m_B} p^{\beta_{m_B}} u(t) & \end{aligned} \quad (6)$$

where $(a_i, b_j)$ are real numbers and $\alpha_1 < \alpha_2 < \dots < \alpha_{m_A}$ and $\beta_0 < \beta_1 < \dots < \beta_{m_B}$ are allowed to be non-integer positive numbers.

The Laplace transform allows to relate the input and the output through a fractional order transfer function model:

$$G(s) = \frac{B(s)}{A(s)} = \frac{\sum_{i=0}^{m_B} b_i s^{\beta_i}}{1 + \sum_{j=1}^{m_A} a_i s^{\alpha_i}}. \quad (7)$$

When all differentiation orders are integer multiple of a basic order $\nu$, called a commensurate order, transfer function (7) becomes

$$G(s) = \frac{B(s)}{A(s)} = \frac{\sum_{i=0}^{m_B} b_i s^{i\nu}}{1 + \sum_{j=1}^{m_A} a_i s^{j\nu}}. \quad (8)$$

Stability of fractional order systems has been analyzed in different contexts. The most well-known stability criterion was established by [49] and allows to check the stability of $\nu$-commensurate order systems 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. Orders with commensurate orders beyond 2 can be proven to be unstable [50]. Further extensions have been developed in order to check stability of non-commensurate systems [51].

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

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

and, for every pole $s_k$ ($P(s_k) = 0$):

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

□

<footnote n="1">¹Note that a too small sampling period may also provide numerical issues such as stability in digital implementation. In this case, suitable discrete rational approximations might be preferred.</footnote>
## B. Problem formulation

Consider an input $u(t)$ and its noise-free output $y(t)$ related by the following compact model form:

$$G(p, \theta) = \frac{B(p)}{A(p)} = \frac{\sum_{i=0}^{m_B} b_i p^{\beta_i}}{1 + \sum_{j=1}^{m_A} a_j p^{\alpha_j}} \quad (11)$$

where $A(p)$ and $B(p)$ are assumed to be coprime and the system stable. Data is collected at samples $T_s, 2T_s, \dots$, for $t = 0$ to the current time $t = kT_s$. The output data is supposed to include an additive white noise $\xi(t)$ due to experiment and sensor imperfections in a realistic scenario. The complete equation can be written in the following form:

$$\begin{cases} y(t) = G(p, \theta)u(t) \\ y^*(t) = y(t) + \xi(t). \end{cases} \quad (12)$$

When using model (11), the parameter vector $\theta$:

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

gathers $\rho (m_a + m_b + 1)$ transfer function coefficients:

$$\rho = [b_0, b_1, \dots, b_{m_B}, a_1, a_2, \dots, a_{m_A}]^T \quad (14)$$

and $\mu (m_a + m_b + 1)$ differentiation orders:

$$\mu = [\beta_0, \beta_1, \dots, \beta_{m_B}, \alpha_1, \alpha_2, \dots, \alpha_{m_A}]^T. \quad (15)$$

The total number of parameters to be estimated is then $2 \times (m_a + m_b + 1)$. If $m_a$ or $m_b$ are too high, 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 (8) can be chosen, where $\mu$ boils down to a single differentiation order $\nu$. As a consequence, the parameter vector $\theta$ is then reduced to $m_a + m_b + 2$ parameters. In the following, only commensurate order model will be studied.

# III. RECURSIVE COEFFICIENT ESTIMATION FOR FRACTIONAL ORDER SYSTEMS

If all differentiation orders are assumed to be known, only the coefficients are to be estimated. Consequently, the parameter vector $\theta$ boils down to:

$$\theta = \rho = [b_0, b_1, \dots, b_{m_B}, a_1, a_2, \dots, a_{m_A}]^T. \quad (16)$$

## A. 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 (17)$$

where the regression vector is as follows:

$$\phi^*(t)^T = \begin{bmatrix} p^{\beta_0} u(t), \dots, p^{\beta_{m_B}} u(t), \\ -p^{\alpha_1} y^*(t), \dots, -p^{\alpha_{m_A}} y^*(t) \end{bmatrix}. \quad (18)$$

An error function can be defined as:

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

By minimizing the quadratic error function:

$$J(\hat{\theta}) = \frac{1}{2} \sum_{h=1}^{k} \epsilon(hT_s, \hat{\theta})^2, \quad (20)$$

the parameter vector $\hat{\theta}$ is its solution according to:

$$\hat{\theta}(t) = \arg \min \sum_{h=1}^{k} \epsilon(hT_s, \hat{\theta})^2. \quad (21)$$
In presence of high level noise, the least-squares leads to biased estimates. In order to reduce the bias, state variable filters (SVF) are used to provide filtered input $u_f(t)$ and filtered output $y_f(t)$. The filtering reduces noise estimation on the estimation and is typically a low-pass filter tuned to a precise bandwidth. Poisson's filter is a common choice:

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

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

$$N_f > \left[ \max(\beta_{m_B}, \alpha_{m_A}) \right], \qquad (23)$$

where $[\cdot]$ is the ceiling operator. The choice of the cut-off frequency $\omega_f$ should not be neglected, as it can have an influence on the convergence rate of the algorithm [52]. Normally, this frequency is chosen to be slightly larger to the system bandwidth [40]. It should be noted that some empirical knowledge of the plant may be required in order to refine it.

The filtered error $\epsilon_f(t)$:

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

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_B}} u_f(t), \\ -p^{\alpha_1} y_f^*(t), \dots, -p^{\alpha_{m_A}} y_f^*(t) \end{bmatrix}. \qquad (25)$$

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} \qquad (26)$$

Inspired by the recursive least-squares [24], the recursive least squares state variable filter (RLSSVF) algorithm is obtained at each iteration $k$ for fractional order systems by following:

$$\left\{ \begin{aligned} 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)} [I - L(k)\phi_f^*(kT_s)^T] F(k-1) \\ \hat{\theta}(k) &= \hat{\theta}(k-1) + L(k) [y_f^*(kT_s) - \phi_f^*(kT_s)^T \hat{\theta}(k-1)] \end{aligned} \right. \qquad (27)$$

where $\gamma$ is a forgetting factor that gives more or less weight to the most recent past of signals and $F$ is the covariance matrix. Remark 3.1: The forgetting factor $\gamma$ may actually be an adaptive function of time according to the iteration [53].

Remark 3.2: In the case of fractional order system, the forgetting factor is useless as the Grünwald-Letnikov fractional derivative (1) owns by itself a forgetting factor: the Newton's binomial coefficient (2). Therefore, in the whole study, the forgetting factor $\gamma$ will be taken as unity.

## B. Recursive Prediction Error Method (RPEM)

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

$$\epsilon(t, \hat{\theta}) = y(t) - \hat{y}(t), \qquad (28)$$

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

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

The following quadratic error function is minimized:

$$J(\hat{\theta}) = \frac{1}{2} \sum_{h=1}^{k} \epsilon(t, \hat{\theta})^2. \qquad (30)$$

By taking the gradient of the criterion:

$$\frac{\partial J(\hat{\theta})}{\partial \hat{\theta}} = - \sum_{h=1}^{k} \psi_{\rho}(hT_s, \hat{\theta}) \epsilon(hT_s, \hat{\theta}) \qquad (31)$$

with

$$\psi_{\rho}(hT_s, \hat{\theta}) = \frac{\partial \hat{y}(hT_s, \hat{\theta})}{\partial \hat{\theta}}, \qquad (32)$$
its error sensitivity is computed by:

$$ \frac{\partial \epsilon(hT_s, \hat{\theta})}{\partial \hat{\theta}} = - \left[ \frac{\partial \hat{y}}{\partial \hat{b}_0}, \dots, \frac{\partial \hat{y}}{\partial \hat{b}_{m_B}}, \frac{\partial \hat{y}}{\partial \hat{a}_1}, \dots, \frac{\partial \hat{y}}{\partial \hat{a}_{m_A}} \right]^T \bigg|_{t=hT_s} \quad (33) $$

where

$$ \frac{\partial \hat{y}}{\partial \hat{b}_i} = \frac{p^{\beta_i}}{1 + \sum_{j=1}^{N} \hat{a}_j p^{\alpha_j}} u(hT_s) \quad (34) $$

and

$$ \frac{\partial \hat{y}}{\partial \hat{a}_{\ell}} = - \frac{\sum_{i=0}^{m_B} \hat{b}_i p^{\beta_i + \alpha_{\ell}}}{\left(1 + \sum_{j=1}^{m_A} \hat{a}_j p^{\alpha_j}\right)^2} u(hT_s). \quad (35) $$

for $i = 0, \dots, m_B$ and $\ell = 1, \dots, m_A$.

Inspired by the recursive least-squares proposed for integer order systems [24], the recursive prediction error method (RPEM) algorithm is obtained for fractional order systems by following these equations at each iteration $k$:

$$ \left\{ \begin{aligned} \epsilon(kT_s, \hat{\theta}(k-1)) &= y^*(kT_s) - \hat{y}(kT_s) \\ R(k) &= R(k-1) + \gamma_{\rho} \left[ \psi_{\rho}(kT_s, \hat{\theta}(k-1)) \right. \\ & \qquad \left. \times \psi_{\rho}^{T}(kT_s, \hat{\theta}(k-1)) - R(k-1) \right] \\ \hat{\theta}(k) &= \hat{\theta}(k-1) + \gamma_{\rho} R^{-1}(k) \\ & \qquad \times \psi_{\rho}(kT_s, \hat{\theta}(k-1)) \epsilon(kT_s, \hat{\theta}(k-1)) \end{aligned} \right. \quad (36) $$

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 [54]. For the present work, it will be taken as a constant.

### C. 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, \hat{\theta}) = [\epsilon(0, \hat{\theta}), \epsilon(T_s, \hat{\theta}), \dots, \epsilon(kT_s, \hat{\theta})]^T. \quad (37) $$

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}(kT_s, \hat{\theta})$ is now a matrix:

$$ \tilde{\psi}_{\rho}(kT_s, \hat{\theta}) = - \begin{bmatrix} \frac{\partial \epsilon(0, \hat{\theta})}{\partial \hat{b}_0} & \frac{\partial \epsilon(T_s, \hat{\theta})}{\partial \hat{b}_0} & \dots & \frac{\partial \epsilon(kT_s, \hat{\theta})}{\partial \hat{b}_0} \\ \frac{\partial \epsilon(0, \hat{\theta})}{\partial \hat{b}_1} & \frac{\partial \epsilon(T_s, \hat{\theta})}{\partial \hat{b}_1} & \dots & \frac{\partial \epsilon(kT_s, \hat{\theta})}{\partial \hat{b}_1} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \epsilon(0, \hat{\theta})}{\partial \hat{a}_{m_A}} & \frac{\partial \epsilon(T_s, \hat{\theta})}{\partial \hat{a}_{m_A}} & \dots & \frac{\partial \epsilon(kT_s, \hat{\theta})}{\partial \hat{a}_{m_A}} \end{bmatrix} \quad (38) $$

where $k$ indicates the present iteration.

Let us introduce the extended measured output $\tilde{Y}^*(kT_s)$ and the extended estimated output $\tilde{Y}(kT_s)$ as:

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

$$ \tilde{Y}(kT_s) = [\hat{y}(0), \hat{y}(T_s), \dots, \hat{y}(kT_s)]^T. \quad (40) $$
The long memory recursive prediction error method (LMRPEM) is proposed for fractional order systems by following these equations at each iteration k:

$$
\left\{
\begin{aligned}
\tilde{\epsilon}(kT_s, \hat{\theta}) &= \tilde{Y}^*(kT_s) - \tilde{Y}^*(kT_s) \\
R(k) &= R(k-1) + \gamma_{\rho} \left[ \tilde{\psi}_{\rho}(kT_s, \hat{\theta}(k-1)) \right. \\
&\qquad \left. \times \tilde{\psi}_{\rho}^T(kT_s, \hat{\theta}(k-1)) - R(k-1) \right] \\
\hat{\theta}(k) &= \hat{\theta}(k-1) + \gamma_{\rho} R^{-1}(k) \\
&\qquad \times \tilde{\psi}_{\rho}(kT_s, \hat{\theta}(k-1)) \tilde{\epsilon}(kT_s, \hat{\theta}(k-1)).
\end{aligned}
\right.
\quad (41)
$$

D. Recursive Simplified Refined Instrumental Variable for Continuous-time Fractional systems (RSRIVCF)

The Simplified Refined Instrumental Variable for Continuous-time (SRIVC) algorithm is derived from SRIV for offline discrete-time model identification. The extension to fractional systems (SRIVCF) was proposed in [31]. A maximum likelihood estimation was evolved to a pseudo-linear form involving optimal pre-filters [55]. Algebraic and statistical similarity allows us to use such concepts in the particular scope of this paper. Following the prediction error method, the error function ϵ can be defined by the output error:

$$
\epsilon(t, \hat{\theta}) = y^*(t) - \frac{B(p, \hat{\theta})}{A(p, \hat{\theta})} u(t) \quad (42)
$$

which can also be rewritten as:

$$
\epsilon(t, \hat{\theta}) = A(p, \hat{\theta}) \left( \frac{1}{A(p, \hat{\theta})} y^*(t) \right) - B(p, \hat{\theta}) \left( \frac{1}{A(p, \hat{\theta})} u(t) \right) \quad (43)
$$

$$
= A(p, \hat{\theta}) y_f^*(t) - B(p, \hat{\theta}) u_f(t) \quad (44)
$$

where

$$
u_f(t) = \frac{1}{A(p, \hat{\theta})} u(t) \quad \text{and} \quad y_f^*(t) = \frac{1}{A(p, \hat{\theta})} y^*(t). \tag{45}
$$

The latter result shows that such a filtering allows going from an output error form (42) to an equation error one (44). The optimal filter that leads to the optimal instrumental variable estimates is given by:

$$
F^{\text{opt}}(p) = \frac{1}{A(p, \theta)} \tag{46}
$$

where θ stands for the true parameter vector. In practical cases, A(p, θ) is unknown, and its estimation is iteratively computed:

$$
F^{\text{iter}}(p) = \frac{1}{A(p, \hat{\theta})} = \frac{1}{1 + \sum_{j=1}^{m_A} \hat{a}_i p^{\alpha_i}}. \quad (47)
$$

The noise-free output estimates are obtained from an auxiliary model whose coefficients are also iteratively computed:

$$
\hat{x}_{IV}^{iter}(t) = \frac{B(p, \hat{\theta})}{A(p, \hat{\theta})} u(t). \quad (48)
$$

An instrumental vector ζ(t, ˆθ) is defined as:

$$
\zeta(t, \hat{\theta}) = \begin{bmatrix} p^{\beta_0} u(t), \dots, p^{\beta_{m_B}} u(t), \\ -p^{\alpha_1} \hat{x}_{IV}^{iter}(t), \dots, -p^{\alpha_{m_A}} \hat{x}_{IV}^{iter}(t) \end{bmatrix}^T . \quad (49)
$$

The filtered instrument ζ_f(t, ˆθ) is then:

$$
\zeta_f(t, \hat{\theta}) = F^{\text{iter}}(p) \times \left[ \begin{array}{l} p^{\beta_0} u(t), \dots, p^{\beta_{m_B}} u(t), \\ -p^{\alpha_1} \hat{x}_{\text{IV}}^{\text{iter}}(t), \dots, -p^{\alpha_{m_A}} \hat{x}_{\text{IV}}^{\text{iter}}(t) \end{array} \right]^T . \quad (50)
$$

The error function is given by the output error:

$$
\epsilon(t, \hat{\theta}) = y(t) - \hat{x}_{IV}^{iter}(t), \quad (51)
$$

and the quadratic error function is given by:

$$
J(\hat{\theta}) = \frac{1}{2} \sum_{h=1}^{k} \epsilon(t, \hat{\theta})^2. \qquad (52)
$$
Finally, the recursive simplified refined instrumental variable for continuous-time fractional systems (RSRIVCF) algorithm is provided for fractional order systems by following these equations at each iteration k:

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

#### IV. RECURSIVE COEFFICIENT AND DIFFERENTIATION ORDER ESTIMATION

The non-trivial task of fully estimating parameter vector as defined in equation (13) will now be explored. The fractional differentiation are now assumed unknown. Note that contrary to the coefficients, the differentiation orders are not expressed in a linear way such as in (7) or (8). Consequently, a suitable error function needs to be defined to help parameter estimation. Contrary to the integer case, recursively estimating the differentiation orders changes the model order at each iteration. For this study, order estimation will be limited to commensurate order systems, such as defined in equation (8). Therefore, a single differentiation order $\nu$ will be estimated.

##### A. Differentiation order estimation

First, assuming that the coefficients are assumed known and that a fractional model $G(p, \nu)$ is $\nu$-commensurate, only the commensurate order $\nu$ is to be estimated. Consequently, the parameter vector $\theta$ boils down to:

$$ \theta = \mu = \nu. \qquad (54) $$

Even for the simplest models, differentiation orders are nonlinear with respect to the transfer function, thus justifying the use of a prediction error method.

First, a suitable output error can be defined as:

$$ \epsilon(t, \mu) = y^*(t) - \hat{y}(t, \mu). \qquad (55) $$

A correction term is then added to (55) in order to recursively estimate the differentiation orders:

$$ \hat{\mu}(k) = \hat{\mu}(k-1) + \gamma_{\mu} R_{\mu}^{-1} \psi_{\mu}(kT_s, \mu(k-1)) \epsilon(kT_s, \hat{\mu}(k-1)) \qquad (56) $$

where:

$$ \psi_{\mu}(kT_s, \hat{\mu}(k-1)) = - \frac{\partial \epsilon(kT_s, \hat{\mu}(k-1))}{\partial \hat{\mu}} = \frac{\partial \hat{y}}{\partial \hat{\mu}} \Big|_{t=kT_s} . \qquad (57) $$

As differentiation orders heavily impact long-memory behavior of fractional models ([56] [48, Chap. 7]), the LMRPEM extended error is proposed:

$$ \tilde{\epsilon}(kT_s, \hat{\mu}) = [\epsilon(0, \hat{\mu}), \epsilon(T_s, \hat{\mu}), \dots, \epsilon(kT_s, \hat{\mu})]^T, \qquad (58) $$

which turns the error sensitivity into a matrix:

$$ \tilde{\psi}_{\mu}(kT_s, \hat{\mu}) = - \begin{bmatrix} \frac{\partial \epsilon(0, \hat{\mu})}{\partial \hat{\beta}_0} & \frac{\partial \epsilon(T_s, \hat{\mu})}{\partial \hat{\beta}_0} & \dots & \frac{\partial \epsilon(kT_s, \hat{\mu})}{\partial \hat{\beta}_0} \\ \frac{\partial \epsilon(0, \hat{\mu})}{\partial \hat{\beta}_1} & \frac{\partial \epsilon(T_s, \hat{\mu})}{\partial \hat{\beta}_1} & \dots & \frac{\partial \epsilon(kT_s, \hat{\mu})}{\partial \hat{\beta}_1} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \epsilon(0, \hat{\mu})}{\partial \hat{\alpha}_{m_A}} & \frac{\partial \epsilon(T_s, \hat{\mu})}{\partial \hat{\alpha}_{m_A}} & \dots & \frac{\partial \epsilon(kT_s, \hat{\mu})}{\partial \hat{\alpha}_{m_A}} \end{bmatrix}. \qquad (59) $$

For commensurate order estimation, the gradient becomes:

$$ \tilde{\psi}_{\mu}(kT_s, \hat{\mu}) = \tilde{\psi}_{\nu}(kT_s, \hat{\nu}) \\ = - \left[ \frac{\partial \epsilon(0, \hat{\nu})}{\partial \hat{\nu}}, \frac{\partial \epsilon(T_s, \hat{\nu})}{\partial \hat{\nu}}, \dots, \frac{\partial \epsilon(kT_s, \hat{\nu})}{\partial \hat{\nu}} \right]. \qquad (60) $$

The error sensitivity function is defined by:

$$ \frac{\partial \epsilon(hT_s, \hat{\nu})}{\partial \nu} = - \frac{\partial \hat{y}}{\partial \nu} \Big|_{t=hT_s}, \qquad (61) $$
**Algorithm 1:** Hybrid LMRPEM-2 algorithm

Step 1:

1) Initialization

k = 0

Initialize parameter vector $\theta(0) = \begin{bmatrix} \rho(0) \\ \nu(0) \end{bmatrix}$

Compute the error function $J(0)$ from (30)

2) Recursive algorithm

for $k = 1 : end$

a) Coefficient estimation by LMRPEM

Assuming the order $\hat{\nu}(k-1)$, compute $\hat{\rho}(k)$ from equation (41)

Compute the error function $J(k)$ from (30)

if $J(k) > J(k-1)$

$$ \gamma_{\rho} = \frac{\gamma_{\rho}}{2} $$

Compute $\hat{\rho}(k)$ from equation (41)

Compute the error function $J(k)$ from (30)

b) Differentiation order estimation by LMRPEM

Assuming the coefficient $\hat{\rho}(k)$ known, compute $\hat{\nu}(k)$ from equation (63)

Compute the error function $J(k)$ from (30)

if $J(k) > J(k-1)$

$$ \gamma_{\nu} = \frac{\gamma_{\nu}}{2} $$

Compute $\hat{\nu}(k)$ from equation (63)

Compute the error function $J(k)$ from (30)

where the output sensitivity $\frac{\partial \hat{y}}{\partial \hat{\nu}}$ is expressed by:

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

Even though this expression remains linear, computing this gradient in the time-domain is non trivial because of the $\ln(p)$ term. In order to simplify calculations, the gradient will be numerically computed instead of analytically. LMRPEM is then used for differentiation order estimation:

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

## B. Coefficient and differentiation order estimation: hybrid algorithms

Let us now consider that both coefficients and differentiation orders are unknown. The differentiation order and coefficient estimation is done by using two-stage algorithms. Two-stage algorithms have been successfully used in the literature in different contexts (see. e.g. [31], [33], [57]–[60]). Although they do not have a proof of convergence, it has been noted in these references that they do often converge to a minimum. They both share the differentiation order estimation stage, as for coefficient estimation, this stage is carried out with two different methods: LMRPEM (see section III-C) and RSRIVC (see section III-D).

In spite of algorithm similarity, differences regarding parameters to be tuned should not be neglected. As stated in [29], the tuning of the forgetting factor $\gamma$ in recursive PEM methods is non-trivial and highly influences the convergence of the parameters. In some cases, methods implying a time-varying factor have been proposed, as seen in [24]. In addition, these classic methods imply the tuning of a single parameter, as orders are not estimated for integer order (or rational) systems. Consequently, algorithm 1 requires the tuning of two forgetting factors (one for the coefficients and one for the differentiation orders), whereas a single one is required for algorithm 2. For this study, the tuning of the parameters was performed empirically by trial and error. If these factors are too significant, convergence towards a minimum may become difficult once the system is near a minimum. Therefore, if an error criterion increment is detected in an iteration, the forgetting factors are reduced and the previous calculated parameters are kept. A difference in the output estimation $\hat{y}(t)$ should also be noted, as the RSRIVCF
**Algorithm 2:** Hybrid LMRPEM-IV algorithm

Step 1:

1) Initialization

$k = 0$

Initialize parameter vector $\theta(0) = \begin{bmatrix} \rho(0) \\ \nu(0) \end{bmatrix}$

Compute the error function $J(0)$ from (30)

2) Recursive algorithm

for $k = 1 : end$

a) Coefficient estimation by RSRIVCF

Assuming the order $\hat{v}(k-1)$ known,
compute $\hat{\rho}(k)$ from equation (53)

Compute the error function $J(k)$ from (52)

b) Differentiation order estimation by LMRPEM

Assuming the coefficient $\hat{\rho}(k)$ known,
compute $\hat{v}(k)$ from equation (63)

Compute the error function $J(k)$ from (30)

if $J(k) > J(k-1)$

$$ \gamma_{\nu} = \frac{\gamma_{\nu}}{2} $$

Compute $\hat{v}(k)$ from equation (63)

Compute the error function $J(k)$ from (30)

algorithm will use the previous measured data whereas the LMRPEM algorithm will continuously apply a convolution product as in equation (29).

## V. SIMULATIONS AND RESULTS

### A. True system description

Let us consider the following true system model:

$$ G(s) = \frac{K}{1 + 2\zeta \left( \frac{s}{\omega_0} \right)^{\nu} + \left( \frac{s}{\omega_0} \right)^{2\nu}} \quad (64) $$

Such a general structure represents a second kind fractional order system [50]. The true parameters are given by $K = 1$, $\zeta = -0.5$, $\omega_0 = 0.5 \text{ rad/s}$ and $\nu = 0.5$.

The general form of model (64) is given by:

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

where true values are $b_0 = K = 1$, $a_2 = \frac{1}{\omega_0^2\nu} = 2$ and $a_1 = \frac{2\zeta}{\omega_0^\nu} = -1.414$. The input is a pseudo-random binary signal oscillating between -5 and 5 and containing 1000 points with a sampling time of $T_s = 0.01 \text{ s}$. According to the system description (12), the output $y(t)$ is corrupted by a zero-mean white Gaussian noise $\xi(t)$ with a signal to noise ratio (SNR) level defined as:

$$ SNR = 10 \log \frac{\|y(t)\|_2^2}{\|\xi(t)\|_2^2} \quad (66) $$

All simulations are initialized with $\rho_0^T = [1, 0.5, 2.5]$ and $\nu_0 = 0.35$ (for hybrid algorithms). The Bode diagrams of the true system and input-output data for simulation are respectively shown in figures 1 and 2.

All simulations have been carried out by using the CRONE toolbox developed in Matlab®, which is dedicated to fractional calculus, fractional system simulation and system identification with fractional models<fnref n="2" />.

### B. Coefficient estimation results

For this subsection, the commensurate order will be supposed as known and taken as equal to its true value $\nu = 0.5$.

From Fig. 1, the true system is resonant in the frequency bandwidth [1, 10]rad/s. Therefore, the Poisson filter (22) will have a cut-off frequency of $\omega_f = 10 \text{ rad/s}$ and as $\max(\beta_{mB}, \alpha_{mA}) = 1$, $N_f = 2$.

<footnote n="2">²The new version of the CRONE toolbox, released in 2022, 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 at http://cronetoolbox.ims-bordeaux.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Bode diagram of the true system (64)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Input/output data used for recursive system identification with SNR = 15dB</figcaption></figure>

Note that in a realistic scenario such information may not be available and an empirical estimation of the bandwidth is required. Otherwise, the input signal should be sufficiently persistently exciting, meaning that its spectrum should have at least half the number of frequency rays than the total number of parameters [61]. The reader can have a look at [62] for experiment planning of fractional systems of the second kind. Without system *a priori* knowledge, the best input signal is a PRBS sequence that has almost a uniform spectral density such as a white noise, but with amplitudes sufficiently high to differentiate it from output measurement noise.

Monte Carlo simulation results for 100 runs are summarized in Table I. As expected, the RLSSVF method is the most sensitive to noise. It has an important estimation standard deviation (and consequently estimation variance) level and estimation bias is already present from the lowest noise level and gets worse as the noise level increases. The algorithm is not efficient as it fails to provide accurate results for the parameter estimation. One way to improve such a poor performance would be to increase to better refine the cut-off frequency $\omega_f$ and take it closer to the system bandwidth. However, it should be noted that such a strategy could deteriorate the estimation in practice as the system bandwidth is usually not available.

The RPEM method ($\gamma_\rho$ being 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 parameter true value is always included in the band determined by the standard
<table><caption>TABLE I: 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></th><th></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>True</th><th></th><th colspan="2">1.000</th><th colspan="2">2.000</th><th colspan="2">-1.414</th></tr><tr><th>Method</th><th>SNR</th><th>ê<sub>0</sub></th><th>σ̂<sub>b<sub>0</sub></sub></th><th>ễ<sub>2</sub></th><th>σ̂<sub>a<sub>2</sub></sub></th><th>ễ<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><tr><td rowspan="3">RSRIVC</td><td>20</td><td>1.010</td><td>0.006</td><td>1.998</td><td>0.007</td><td>-1.392</td><td>0.008</td></tr><tr><td>15</td><td>1.010</td><td>0.009</td><td>1.999</td><td>0.012</td><td>-1.393</td><td>0.013</td></tr><tr><td>10</td><td>1.008</td><td>0.017</td><td>2.003</td><td>0.025</td><td>-1.393</td><td>0.026</td></tr></tbody></table>

deviation. On the other hand, a rough factor of 3 can be observed between the estimation standard deviations for *SNR* = 20dB and *SNR* = 10dB. Such as the RLSSVF method, the deviations are amongst the most significant ones.

The LMRPEM method ($\gamma_{\rho}$ being set to 0.01) has the most consistent estimation as the bias is negligible 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. In mean, the estimated parameter are very close to the true parameters.

Finally, the RSRIVC method proves to be in an intermediate position between the LMRPEM and RPEM methods. As can be observed, this method provides similar estimation standard deviations to the ones obtained with the LMRPEM method. A slight increase may be observed when the noise level increases (when the *SNR* decreases), but the method still manages to significantly reduce the estimation standard deviation with respect to RLSSVF and RPEM methods. The true parameter values are still contained within the limits of the standard deviations, but the estimation of the parameter $a_1$ is significantly closer to its lower limit when compared to the results with LMRPEM method. The main advantage with this method is the lack of a parameter $\gamma_{\rho}$.

### C. All parameter estimation: differentiation order and coefficient estimation with hybrid algorithms

For this section, the tuning of the $\gamma_X$ gains was performed empirically and is given in Table II.

<table><caption>TABLE II: Hybrid algorithm parameters</caption><tr><td>LMRPEM-2 algorithm</td><td>$\gamma_{\rho} = 0.005$ and $\gamma_{\nu} = 0.05$</td></tr><tr><td>LMRPEM-IV algorithm</td><td>$\gamma_{\nu} = 0.05$</td></tr></table>

The parameter evolution for an *SNR* = 20dB are respectively plotted in Fig. 3 for hybrid LMRPEM-2 algorithm (see algorithm 1) and Fig. 4 for LMRPEM-IV algorithm (see algorithm 2). Convergence rates are similar for both algorithms. Slight oscillations are visible in the parameter estimation with both hybrid algorithms during the transitory phase. Once the differentiation order estimation has converged, the coefficient estimation gets better. Note that both algorithms tend to the commensurate order true value of $\nu = 0.5$.

<table><caption>TABLE III: Monte Carlo simulation results (100 runs) for all parameter estimation with different *SNR*. $\tilde{\theta}$ denotes the parameter mean value and $\hat{\sigma}_{\theta}$ the standard deviations</caption><thead><tr><th></th><th></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><th colspan="2">v</th></tr><tr><th>True</th><th></th><th colspan="2">1.000</th><th colspan="2">2.000</th><th colspan="2">-1.414</th><th colspan="2">0.5</th></tr><tr><th>Method</th><th>SNR</th><th>ê<sub>0</sub></th><th>σ̂<sub>b<sub>0</sub></sub></th><th>ễ<sub>2</sub></th><th>σ̂<sub>a<sub>2</sub></sub></th><th>ễ<sub>1</sub></th><th>σ̂<sub>a<sub>1</sub></sub></th><th>ṽ</th><th>σ̂<sub>v</sub></th></tr></thead><tbody><tr><td>Hybrid</td><td>∞</td><td>0.991</td><td></td><td>2.001</td><td></td><td>-1.431</td><td></td><td>0.497</td><td></td></tr><tr><td rowspan="2">LMRPEM-2</td><td>20</td><td>0.987</td><td>0.023</td><td>2.001</td><td>0.007</td><td>-1.437</td><td>0.042</td><td>0.496</td><td>0.007</td></tr><tr><td>15</td><td>0.989</td><td>0.029</td><td>2.001</td><td>0.010</td><td>-1.434</td><td>0.052</td><td>0.496</td><td>0.009</td></tr><tr><td>Hybrid</td><td>∞</td><td>0.992</td><td></td><td>2.013</td><td></td><td>-1.430</td><td></td><td>0.499</td><td></td></tr><tr><td rowspan="2">LMRPEM-IV</td><td>20</td><td>0.956</td><td>0.045</td><td>2.015</td><td>0.008</td><td>-1.497</td><td>0.084</td><td>0.487</td><td>0.015</td></tr><tr><td>15</td><td>0.934</td><td>0.085</td><td>2.012</td><td>0.015</td><td>-1.531</td><td>0.153</td><td>0.480</td><td>0.028</td></tr></tbody></table>

Monte Carlo Simulation are provided in Table III. Three levels of noise are considered and 100 runs are carried out.

In the noiseless scenario ($SNR = \infty$), both hybrid algorithms provide similar results. As the noise level increases (*SNR* decreases), the LMRPEM-2 presents more consistent estimates, with minimum estimation standard deviations (and consequently
<figure><img src="image_4.png" /><figcaption>Fig. 3: Parameter estimation evolution with hybrid LMRPEM-2 (see algorithm 1 for SNR = 20dB</figcaption></figure>

estimation variances) and without bias. Note that, as the commensurate order has been well estimated (without bias), the coefficients are consequently well estimated. With the LMRPEM-IV, the standard deviations are higher than with the previous method, and mean parameter values present a bias as the noise level increases (as SNR decreases). Consequently, the LMRPEM-IV is more sensitive to noise. It can particularly be seen that, as the commensurate order estimation is biased, the coefficient estimation tries to compensate this error.

### D. All parameter estimation: effect of varying initialization

The performances of the two-stage algorithms are now studied with different initialization of the commensurate order $v_0$. By choosing a fixed $SNR = 20dB$, 100 Monte Carlo runs are performed on both algorithms in order to see the resulting estimated parameters. Two different initial conditions are chosen: $v_0 = 0.25$ and $v_0 = 0.60$. Coefficient initial values are kept at the same values as in the rest of this section. Monte Carlo Simulation analysis is provided in Table IV.

As can be observed from the results, the initial condition on the commensurate order does not have a high impact on the standard deviations, as they are fairly similar for both cases. However, the mean values exhibit some sensitivity to this initialization, most notably parameter $a_1$. In both cases, the final estimated value is closer to its true value when the initial value is closer to it (in this case $v_0 = 0.60$). LMRPEM-IV algorithm is considerably more sensitive to the initial condition variations compared to LMRPEM-2 one: indeed LMRPEM-IV exhibits significantly higher standard deviations, and its estimated parameters are farther from the true values when compared to the results provided by LMRPEM-2. Consequently, the LMRPEM-2 algorithm seems more robust to the initial condition variations.

<table><caption>TABLE IV: Monte Carlo simulation results (100 runs) for all parameter estimation with different $v_0$. $\bar{\theta}$ denotes the parameter mean value and $\hat{\sigma}_{\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><td>True</td><td>1.000</td><td>2.000</td><td>-1.414</td><td>0.5</td></tr><tr><td>Method</td><td>ν<sub>0</sub></td><td>β̂<sub>0</sub></td><td>σ̂<sub>b<sub>0</sub></sub></td><td>σ̂<sub>a<sub>2</sub></sub></td><td>σ̂<sub>a<sub>1</sub></sub></td><td>ν̂</td><td>σ̂<sub>ν</sub></td></tr><tr><td>Hybrid</td><td>0.25</td><td>0.982</td><td>0.017</td><td>2.001</td><td>0.006</td><td>-1.447</td><td>0.032</td><td>0.494</td><td>0.005</td></tr><tr><td>LMRPEM-2</td><td>0.60</td><td>0.997</td><td>0.016</td><td>1.999</td><td>0.007</td><td>-1.419</td><td>0.029</td><td>0.499</td><td>0.005</td></tr><tr><td>Hybrid</td><td>0.25</td><td>0.887</td><td>0.044</td><td>2.014</td><td>0.008</td><td>-1.620</td><td>0.077</td><td>0.464</td><td>0.015</td></tr><tr><td>LMRPEM-IV</td><td>0.60</td><td>1.072</td><td>0.044</td><td>1.991</td><td>0.021</td><td>-1.268</td><td>0.092</td><td>0.524</td><td>0.014</td></tr></tbody></table>

## VI. APPLICATION: THERMAL IMPEDANCE IN LUNGS

Mechanics in lung and heat transfers may be critical in some scenarios. It is usual for cardiac surgery to require the use of extracorporeal circulation (ECC), which could damage lungs [63]. The use of hypothermia has been proposed to protect lungs [64]. Therefore, lung thermal impedance may be relevant to improve thermal control in such cases.
<figure><img src="image_5.png" /><figcaption>Fig. 4: Parameter estimation evolution with hybrid LMRPEM-IV (see algorithm 2) for SNR = 20dB</figcaption></figure>

## A. Lung model

A global thermal model is proposed in [65] by using a series of thermal two-port networks (Fig. 5).

Each single section of the lung is modelled by using a thermal two-port network [10], as shown in Fig. 6. Lung blood flow is non-uniform and the first branches are considered as a dead space used for passive air transport. Blood perfusion becomes important around branch 15-th. Therefore, a mixed model will be used: the first 15 ($n < 15$) branches will be considered as air, whereas the rest of the lung will be treated as having blood perfusion. The thermo-physical parameters of blood and tissue for a human body can be found in [66] and are presented in Table V, where the subscript $t$ stands for tissue, $b$ for blood and $a$ for air.

<table><caption>TABLE V: Simulation parameters for lung impedance</caption><tr><td>a<sub>a</sub></td><td>2.2e-5 [m<sup>2</sup>][s<sup>-1</sup>]</td><td>k<sub>a</sub></td><td>0.025 [W][m<sup>-1</sup>][K<sup>-1</sup>]</td></tr><tr><td>ρ<sub>t</sub></td><td>550 [kg][m<sup>-3</sup>]</td><td>ρ<sub>b</sub></td><td>1060 [kg][m<sup>-3</sup>]</td></tr><tr><td>c<sub>t</sub></td><td>3718 [J][kg<sup>-1</sup>][K<sup>-1</sup>]</td><td>c<sub>b</sub></td><td>3900 [J][kg<sup>-1</sup>][K<sup>-1</sup>]</td></tr><tr><td>k<sub>t</sub></td><td>0.28 [W][m<sup>-1</sup>][K<sup>-1</sup>]</td><td>ω<sub>b</sub></td><td>0.0043 s<sup>-1</sup></td></tr><tr><td>L<sub>0</sub></td><td>0.1 m</td><td>d<sub>0</sub></td><td>1.4 cm</td></tr></table>

## B. Lung equivalent impedance and transfer function

This section deals with the thermal equivalent impedance of the lung seen from its entry, namely $n = 0$. Fig. 7 shows the frequency response of the equivalent lung impedance up to level 15. A simple first species fractional order model is used to approximate the lung impedance model [67]:

$$Z_{\text{lung},0-15}(s) = \frac{K}{1 + \left(\frac{s}{\omega_c}\right)^{\nu}}, \qquad (67)$$

whose general form is:

$$G(s) = \frac{b_0}{a_1 s^{\nu} + 1} = \frac{70155}{51.23 s^{0.5} + 1} \qquad (68)$$

with the true values $b_0 = 70155$, $a_1 = 51.23$ and $\nu = 0.5$.

## C. Lung equivalent impedance identification

A pseudo random binary sequence is used as input with an amplitude going from -5 to 5mW. The total time length of the signal is of 4000 seconds with a sampling time of $T_s = 2$ s. The output is the temperature variation at $n = 0$ measured in Kelvin
<figure><img src="image_6.png" /><figcaption>Fig. 5: Lung scheme</figcaption></figure>

(K). The first 350 seconds of the input-output data are shown in Fig. 8. An SNR = 20 dB is used for this simulated data. The estimated coefficients are initialized at $\hat{b}_0 = 30000$, $\hat{a}_1 = 10$ and $\nu_0 = 0.5$. The coefficients for LMRPEM-2 algorithm are set to $\gamma_\rho = 0.05$ and $\gamma_\nu = 0.05$.

Applying the LMRPEM-2 (algorithm 1) leads to parameter estimation evolution shown in Fig. 9. As can be seen, the estimation converges towards its true values. Indeed, the parameter evolution is slow given the true system form which is non-resonant. At the end of the simulation, the estimated lung model has converged to:

$$ \hat{Z}_{\text{lung}}(s) = \frac{71130}{51.77s^{0.497} + 1}, \qquad (69) $$

which is very close to the true parameters.

The relative errors for parameter $b_0$, $a_1$ and $\nu$ are respectively of -1.39%, -1.04% and 0.60%. It can also be seen in Fig. 9 that the parameters have converged around 1000 iterations which corresponds to 2000s (as the sampling time is $T_s = 2s$), knowing that the system bandwidth is around $3e^{-4}$, leading to a time response around 10 000s. As a consequence, the algorithm LMRPEM-2 has converged in one fifth of the system time response.

## VII. CONCLUSIONS AND PERSPECTIVES

Three direct extensions for recursive system identification of continuous-time fractional order models were proposed for coefficient estimations: the least squares (RLSSVF), prediction error method (RPEM) and instrumental variable (RSRIVCF). An improved RPEM method was introduced to take into account the long-memory behavior of fractional order systems (LMRPEM). 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. However, it was observed that LMRPEM and RSRIVCF provide less estimation standard deviations when compared to RPEM.

Then, for all parameter estimation, especially for differentiation order estimation, two variant algorithms were introduced: the LMRPEM-2 algorithm for both coefficient and order estimation with LMRPEM and the LMRPEM-IV algorithm where RSRIVCF is used for coefficient estimation and the LMRPEM for differentiation order estimation.

The LMRPEM-2 algorithm proved to be more consistent than the LMRPEM-IV as the estimations are without bias and with lowest variances. For RPEM methods, the forgetting factor gain $\gamma_X$ is to be tuned to help the convergence. Also, the longer the dataset, the better the convergence to the true parameters. The lung application has proved the efficiency of the LMRPEM-2 algorithm.
<figure><img src="image_7.png" /><figcaption>Fig. 6: Global circuit model of a lung</figcaption></figure>

Research perspectives may include further studies regarding the tuning of 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 signals.

## REFERENCES

[1] B. Riemann, "Versuch einer allgemeinen Auffassung der Integration und Differentiation," *Gesammelte Mathematische Werke und Wissenschaftlicher*, pp. 331–344, Leipzig 1876.

[2] J. Liouville, "Mémoire sur quelques questions de géométrie et de mécanique et sur un nouveau genre de calcul pour résoudre ces équations," *Journal de l'Ecole Polytechnique*, vol. 13, pp. 71–162, 1832.

[3] K. Oldham and J. Spanier, "The replacement of Fick's laws by a formulation involving semidifferentiation," *J. of Electroanalytical Chemistry and Interfacial Electrochemistry*, vol. 26, no. 2-3, pp. 331–341, 1970.

[4] —, *The fractional calculus - Theory and Applications of Differentiation and Integration to Arbitrary Order*. Academic Press, New-York and London, 1974.

[5] A. Oustaloup, *La dérivation non-entière : théorie, synthèse et applications*. Paris: Hermès, 1995.

[6] J. C. Trigeassou and N. Maamri, "Initial conditions and initialization of linear fractional differential equations," *Signal Processing*, vol. 91, no. 3, pp. 427–436, 2011.

[7] E. R. Weibel, "Mandelbrot's fractals and the geometry of life: A tribute to benoît mandelbrot on his 80th birthday," in *Fractals in Biology and Medicine*, G. A. Losa, D. Merlini, T. F. Nonnenmacher, and E. R. Weibel, Eds. Basel: Birkhäuser Basel, 2005, pp. 3–16.

[8] S. Havriliak and S. Negami, "A complex plane representation of dielectric and mechanical relaxation processes in some polymers," *Polymer*, vol. 8, pp. 161–210, 1967.

[9] O. Cois, A. Oustaloup, T. Poinot, and J.-L. Battaglia, "Fractional state variable filter for system identification by fractional model," in *2001 European Control Conf. (ECC)*, Porto, Portugal, 2001, pp. 2481–2486.

[10] J. Duhe, S. Victor, P. Melchior, Y. Abdelmoumen, and F. Roubertie, "Modeling thermal systems with fractional models: human bronchus application," *Nonlinear Dynamics*, vol. 108, pp. 579–595, 2022.

[11] H. Woodson and J. Melcher, *Electromechanical Dynamics-Part I: Discrete Systems*, 05 2012.

[12] B. Krishna, "Studies on fractional order differentiators and integrators: A survey," *Signal Processing*, vol. 91, no. 3, pp. 386–426, 2011.

[13] M. Nakagawa and K. Sorimachi, "Basic characteristics of a fractance device," *IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences*, vol. 75, pp. 1814–1819, 1992.

[14] F. Yakuphanoglu, B. Krishna, and K. Reddy, "Active and passive realization of fractance device of order 1/2," *Active and Passive Electronic Components*, vol. 2008, p. 369421, 2008.

[15] A. Elwakil, "Fractional-order circuits and systems: An emerging interdisciplinary research area," *IEEE Circuits and Systems Magazine*, vol. 10, no. 4, pp. 40–50, 2010.

[16] R. Magin and M. Ovadia, "Modeling the cardiac tissue electrode interface using fractional calculus," *2nd IFAC Workshop on Fractional Differentiation and its Applications*, vol. 39, no. 11, pp. 302–307, 2006.

[17] C. Ionescu, D. Copot, and R. De Keyser, "Respiratory impedance model with lumped fractional order diffusion compartment," *IFAC Proceedings Volumes*, vol. 46, no. 1, pp. 260–265, 2013, 6th IFAC Workshop on Fractional Differentiation and Its Applications.

[18] D. Copot, R. De Keyser, E. Derom, and C. Ionescu, "Structural changes in the copd lung and related heterogeneity," *PLOS ONE*, vol. 12, no. 5, 05 2017.

[19] M. Yu, Y. Li, I. Podlubny, F. Gong, Y. Sun, Q. Zhang, Y. Shang, B. Duan, and C. Zhang, "Fractional-order modeling of lithium-ion batteries using additive noise assisted modeling and correlative information criterion," *J. of Advanced Research*, vol. 25, pp. 49–56, 2020.

[20] A. Ali, A. Radwan, and A. Soliman, "Fractional order butterworth filter: Active and passive realizations," *IEEE J. on Emerging and Selected Topics in Circuits and Systems*, vol. 3, no. 3, pp. 346–354, 2013.
<figure><img src="image_8.png" /><figcaption>Fig. 7: Gain of $Z_{lung,0-15}(j\omega)$</figcaption></figure>

[21] J. Receveur, S. Victor, and P. Melchior, “New interpretation of fractional potential fields for robust path planning,” *Fractional Calculus & Applied Analysis*, vol. 22, no. 1, pp. 113–127, 2019.

[22] —, “Autonomous car decision making and trajectory tracking based on genetic algorithms and fractional potential fields,” *Intelligent Service Robotics*, vol. 13, pp. 315–330, 02 2020.

[23] S. Victor, J. Duhé, P. Melchior, Y. Abdelmoumen, and F. Roubertie, “Long memory recursive prediction error method for identification of continuous-time fractional models,” *Nonlinear Dynamics*, vol. 110, pp. 635–648, 06 2022.

[24] L. Ljung, *System identification – Theory for the user*, 2nd ed. Upper Saddle River, N.J., USA: Prentice-Hall, 1999.

[25] C. Zeng and S. Liang, “Comparative study of discretization zero dynamics behaviors in two multirate cases,” *International Journal of Control, Automation and Systems*, vol. 13, no. 4, pp. 831–842, 2015.

[26] H. Garnier and L. Wang, *Identification of continuous-time models from sampled data*. Springer-Verlag, 2008.

[27] L. Le Lay, “Identification temporelle et fréquentielle par modèle non entier,” Ph.D. dissertation, Université Bordeaux 1, France, 1998.

[28] B. Mathieu, A. Oustaloup, and F. Levron, “Transfer function parameter estimation by interpolation in the frequency domain,” in *European Control Conference, ECC'95*, Rome, Italie, 1995, pp. 2131–2136.

[29] P. Young, “Parameter estimation for continuous-time models – a survey,” *Automatica*, vol. 17, no. 1, pp. 23–29, 1981.

[30] R. Malti, S. Victor, A. Oustaloup, and H. Garnier, “An optimal instrumental variable method for continuous-time fractional model identification,” in *The 17th IFAC World Congress (IFAC'08)*, Seoul, Korea, 2008, pp. 14 379–14 384.

[31] S. Victor, R. Malti, H. Garnier, and A. Oustaloup, “Parameter and differentiation order estimation in fractional models,” *Automatica*, vol. 49, no. 4, pp. 926–935, 2013.

[32] A. Mayoufi, S. Victor, R. Malti, M. Chetouï, and M. Aoun, “Output error MISO system identification using fractional models,” *Fractional Calculus & Applied Analysis*, vol. 5, no. 24, pp. 1601–1618, 2021.

[33] S. Victor, K. Ruiz, P. Melchior, and S. Chaumette, “Dynamical repulsive fractional potential fields in 3d environment,” *Fractional Calculus & Applied Analysis*, pp. 321–345, 04 2022.

[34] R. Malti, J. Sabatier, and H. Akçay, “Thermal modeling and identification of an aluminium rod using fractional calculus,” in *15th IFAC Symposium on System Identification (SYSID'2009)*, St Malo, France, 2009, pp. 958–963.

[35] S. Victor, P. Melchior, R. Malti, and A. Oustaloup, “Robust motion planning for a heat rod process,” *Journal of Nonlinear Dynamics*, vol. 86, no. 2, pp. 1271–1283, 2016.

[36] A. Eddine, B. Huard, J.-D. Gabano, and T. Poinot, “Initialization of a fractional order identification algorithm applied for lithium-ion battery modeling in time domain,” *Communications in Nonlinear Science and Numerical Simulation*, vol. 59, pp. 375–386, 2018.

[37] L. Sommacal, P. Melchior, A. Oustaloup, J.-M. Cabelguen, and A. Ijspeert, “Fractional multi-models of the frog gastrocnemius muscle,” *J. of Vibration and Control*, vol. 14, no. 9-10, pp. 1415–1430, 2008.

[38] A. Oustaloup, F. Levron, S. Victor, and L. Dugard, “Non-integer (or fractional) power model to represent the complexity of a viral spreading: Application to the COVID-19,” *Annual Reviews in Control*, vol. 52, pp. 523–542, 2021.

[39] A. Djouambi, A. V. Besançon, and A. Charef, “Fractional system identification using recursive algorithms approach,” in *2007 European Control Conference (ECC)*, 2007, pp. 1436–1441.

[40] A. Padilla, “Recursive identification of continuous-time systems with time-varying parameters,” Theses, Université de Lorraine, 07 2017.

[41] D. Idiou, A. Charef, A. Djouambi, and A. Voda, “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, November 2014.

[42] S. Victor, P. Melchior, M. Pellet, and A. Oustaloup, “Lung thermal transfer system identification with fractional models,” *IEEE Transactions on Control Systems Technology*, vol. 28, no. 1, pp. 172–182, 2020.

[43] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo, *Theory and Applications of Fractional Differential Equations, Volume 204 (North-Holland Mathematics Studies)*. New York, NY, USA, 2006.

[44] R. Garrappa, E. Kaslik, and M. Popolizio, “Evaluation of fractional integrals and derivatives of elementary functions: Overview and tutorial,” *Mathematics*, vol. 7, no. 407, 2019.

[45] S. Samko, A. Kilbas, and O. Marichev, *Fractional integrals and derivatives: theory and applications*. Gordon and Breach Science, 1993.
<figure><img src="image_9.png" /><figcaption>Fig. 8: Input/output data used for lung equivalent impedance recursive identification (SNR = 20dB)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Evolutions of the lung equivalent impedance parameter with LMRPEM-2 algorithm</figcaption></figure>

<footnote>[46] A. Grünwald, “Über begrenzte Derivationen und deren Anwendung,” *Zeitschrift für Mathematik und Physik*, pp. 441–480, 1867.</footnote>

<footnote>[47] A. Letnikov, “Theory of differentiation of arbitrary order (russian),” *Matematiceskij Sbornik* (Moscou), vol. 3, no. 1, pp. 1–68, 1868.</footnote>

<footnote>[48] I. Podlubny, *Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications*. San Diego: Academic Press, 1999.</footnote>

<footnote>[49] D. Matignon, “Stability properties for generalized fractional differential systems,” *ESAIM proceedings – Systèmes Différentiels Fractionnaires – Modèles, Méthodes et Applications*, vol. 5, 1998.</footnote>

<footnote>[50] R. Malti, X. Moreau, F. Khemane, and A. Oustaloup, “Stability and resonance conditions of elementary fractional transfer functions,” *Automatica*, vol. 47, no. 11, pp. 2462–2467, 2011.</footnote>

<footnote>[51] M. Rivero, S. Rogosin, J. Tenreiro Machado, and J. Trujillo, “Stability of fractional order systems,” *Mathematical Problems in Engineering*, vol. 2013, p. 356215, 2013.</footnote>

<footnote>[52] C. De Wit, “Recursive estimation of the continuous-time process parameters,” in *1986 25th IEEE Conference on Decision and Control*, 1986, pp. 2016–2020.</footnote>

<footnote>[53] X. Sun, J. Ji, B. Ren, C. Xie, and D. Yan, “Adaptive forgetting factor recursive least square algorithm for online identification of equivalent circuit model parameters of a lithium-ion battery,” *Energies*, vol. 12, no. 12, 2019.</footnote>

<footnote>[54] L. Ljung, “Analysis of a general recursive prediction error identification algorithm,” *Automatica*, vol. 17, no. 1, pp. 89–99, 1981.</footnote>
[55] P. Young, Recursive estimation and Time-Series Analysis, 2nd ed. Berlin: Springer-Verlag, 2011.

[56] F. Mainardi, A. Mura, G. Pagnini, and R. Gorenflo, "Fractional relaxation and time-fractional diffusion of distributed order," in 2nd IFAC Workshop on Fractional Differentiation and its Applications, vol. 39, no. 11, 2006, pp. 1–21.

[57] E. Bai, "An optimal two stage identification algorithm for Hammerstein-Wiener nonlin. systems," Automatica, vol. 34, no. 3, pp. 333–338, 1998.

[58] S. Mukhopahay, A. Patra, and G. P. Rao, "Irreducible model estimation for MIMO systems," International Journal of Control, vol. 53, no. 1, pp. 223–253, 1991.

[59] G. Rao and H. Unbehauen, "Identification of continuous-time systems," Electrical Engineering, vol. 153, pp. 185 – 220, 2006.

[60] W. Spinelli, L. Piroddi, and M. Lovera, "A two-stage algorithm for structure identification of polynomial NARX models," American Control Conference Vols 112, vol. 1-12, pp. 2387–2392, 2006.

[61] G. Goodwin and R. Payne, Dynamic system identification. Experiment design and data analysis. Academic Press, 1977.

[62] R. Malti, A. Mayoufi, and S. Victor, "Experiment design for elementary fractional models," Communications in Nonlinear Science and Numerical Simulation, 2022.

[63] R. Badenes, A. Lozano, and F. Belda, "Postoperative pulmonary dysfunction and mechanical ventilation in cardiac surgery," Critical Care Research and Practice, p. 420513, 2015.

[64] R. Giesinger and P. McNamara, "The Impact of Therapeutic Hypothermia on Pulmonary Hemodynamics of Meconium Aspiration Syndrome," Amer. J. of Respiratory and Critical Care Medicine, vol. 198, no. 2, 03 2018.

[65] J.-F. Duhe, S. Victor, P. Melchior, Y. Abdelmoumen, and F. Roubertie, "Global thermal modeling of lung heat transfer with blood perfusion," in International Conference on Fractional Differentiation and its Applications ICFDA 2021, Springer, Ed., vol. 452, Warsaw, Poland, 2021.

[66] A. Coccarelli, E. Boileau, D. Parthimos, and P. Nithiarasu, "An advanced computational bioheat transfer model for a human body with an embedded systemic circulation," Biomechanics and Modeling in Mechanobiology, vol. 15, no. 5, pp. 1173–1190, 2016.

[67] J.-F. Duhe, S. Victor, P. Melchior, Y. Abdelmoumen, and F. Roubertie, "Online system identification of heat transfers in lungs with the LMRPEM-2 method," in 2023 International Conference on Fractional Differentiation and Its Applications (ICFDA), Ajman, United Arab Emirates, 03 2023, pp. 1–6.

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

Jean-François Duhé was born en 1993 in Panama City, Panama. He received his Bachelor degree in electromechanical engineering from the Universidad Tecnologica de Panama (UTP), Panama, in 2016. He completed his MSc in automatic control, from the Université de Bordeaux in 2019 and his PhD in 2022. He is currently an associate professor at Universidad de Panamá. His research interest lies in the area of dynamic systems, fractional calculus, system identification and thermal systems.

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

Stéphane Victor was born in Germany in 1983. He graduated from ENSEIRB and Ecole Polytechnique de Montréal engineering schools in 2006, and obtained his MSc and PhD degrees in Automatic control from Université de Bordeaux in 2006 and 2010. He received his Accreditation to supervise research (HDR) from Univ. Bordeaux in 2022. He is currently Associate Professor and has joined the CRONE team of IMS laboratory at Université de Bordeaux in 2006. His research interests are in fractional calculus, automatic control, system identification, thermal systems, autonomous vehicle, trajectory planning and flatness.

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

Pierre Melchior was born in France on 1961. He graduated from the Univ. Bordeaux, Talence, France, and received the Ph.D. degree in control eng. from Paul Sabatier University, Toulouse, France, in 1989 and the Accreditation to supervise research (HDR) from the Univ. Bordeaux, in 2014. He is currently Full Professor of control theory and researcher in Laboratory IMS, Univ. Bordeaux, Talence. He is head of the CRONE team of IMS Automatic Control Group. His research interests include control theory, robotics, and biomechatronics (e.g., fractional systems, path planning, motion planning, and muscle modeling).
<figure><img src="image_14.png" /></figure>

**Youssef Abdelmoumen** was born in 1967. He obtained his diploma of doctor of general medicine in 1997 and his specialisation in anaesthesia and intensive care domains. He was former hospital practitioner and works at clinique St Augustin, Bordeaux. His research interests are in the area of lungs thermal transfers.

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

**François Roubertie** was born in 1973. Graduated at the medical school of Bordeaux, France in 2002, he received the Ph.D. degree in cell biology and physiopathology from University of Bordeaux, in 2015 and the Accreditation to supervise research (HDR) from the Univ. Bordeaux, in 2016. He is currently Professor and associate director of the Pediatric and Adult Congenital Cardiac Surgery Unit, in Bordeaux Heart Hospital University. He is principal and associate investigator for several clinical and experimental research projects (in the field of congenital heart disease) in both Bordeaux Heart Hospital and in IHU LIRYC - Electrophysiology and Heart Modeling Institute-Bordeaux.