# Online system identification of heat transfers in lungs with the LMRPEM-2 method

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 Roubertie. Online system identification of heat transfers in lungs with the LMRPEM-2 method. 2023 International Conference on Fractional Differentiation and Its Applications (ICFDA), Mar 2023, Ajman, United Arab Emirates. pp.1-6, 10.1109/ICFDA58234.2023.10153153. hal-04317759

HAL Id: hal-04317759

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

Submitted on 12 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Online system identification of heat transfers in lungs with the LMRPEM-2 method

Jean-François Duhé Univ. Bordeaux, CNRS, Bordeaux INP,
IMS, UMR 5218

F-33400 Talence, France Stéphane Victor Univ. Bordeaux, CNRS, Bordeaux INP,
IMS, UMR 5218

F-33400 Talence, France

stephane.victor@ims-bordeaux.fr Pierre Melchior Univ. Bordeaux, CNRS, Bordeaux INP,
IMS, UMR 5218

F-33400 Talence, France Youssef Abdelmoumen Clinique Saint Augustin,
Centre de Recherche Cardio-Thoracique
de Bordeaux (CRCTB) - INSERM U1045

François Roubertie IHU Liryc, Electrophysiology and Heart Modeling Institute,
Fondation Bordeaux Université, 33000, Bordeaux Cedex, France
Centre de Recherche Cardio-Thoracique
de Bordeaux (CRCTB) - INSERM U1045

## Abstract

During cardiac surgery, an extracorporeal circulation system is used. This implies that the lungs are disconnected from the body during the procedure. In order to reduce potential tissue damage, a mild hypothermia is applied to the lungs. An improved comprehension of lung heat transfer dynamics may enhance the control of such hypothermia.

A thermal model based upon the thermal two-port network allows to have a broad frequency domain validity, which allow to have a valid model over the human breathing frequency range. This formalism can be adapted in order to also include the effect of blood flow, which is a natural thermal regulator for the human body. This is achieved through the combination of the thermal two-port network and the bio-heat equation.

Besides the physical modeling, system identification may also be important. If one wants to track model parameter fluctuations during a procedure, a real-time system identification method is desired. Therefore, recursive methods will be developed.

## Index Terms

heat transfer, lungs, heat equation, bio-heat equation, system identification, fractional derivative, fractional calculus

## I. INTRODUCTION

There is a wide range of applications for which temperature may be critical, such as electronic device cooling, electrical machine temperature control or biomedical scenarios. In a general way, heat equation may be applied for any of these cases in order to analyze thermal phenomena (normally through finite difference or finite element methods).

In biology and medicine, temperature fluctuation range is usually very narrow. Even moderate temperature fluctuations may damage organic tissue. In normal conditions, the interior of a human body remains globally close to 37°C. However, cardiac surgery procedures implies organ exposition to the exterior surroundings. As a consequence, heat losses occur. It is indeed one of the main interest of surgeons regards lungs, as they are disconnected from the rest of the body during extracorporeal circulation (ECC). They are exposed to harsh conditions and damage on lung tissue may be important [1]. Mild hypothermia is often employed [2] to reduce metabolism and oxygen consumption, but it is often achieved through empirical methods. Dynamic models of heat transfers may allow an improved temperature control.

Through the last century, an important amount of studies have analyzed the lungs, mostly regarding its elastic properties. Hilderbrandt proposed some early models analyzed its elasticity by performing experiments on balloons as well as cat lungs [3], [4]. The mathematical relationships used by Hildebrandt implied an almost constant-phase behavior in the frequency domain. Even though he didn't directly referred to fractional calculus, the characteristics of these first models indirectly suggest a fractional order dynamic system. Besides this model, simple pneumatic circuit models have also been developed by decades [5]–[7]. They provide a simplified equivalent mechanical impedance of the lungs as a whole system.

Another characteristic regarding lung modeling has also been studied for decades: its geometry. The arborescent structure and recursivity of the lungs was first studied by Mandelbrot [8]. Ionescu and Oustaloup proposed a circuit model which takes into account the particular geometry of the lungs [9], [10]. Even though the resulting model is theoretically interesting, such a model is difficult to use for experimental scenarios. As a consequence, Ionescu proposed 4 fractional order impedance model in order to analyze experimental data. This equivalent impedance models use constant phase elements, which have already been
used for different applications [11], [12]. Once again, the presence of fractional order elements proves to be useful. Fluctuations on fractional order impedance parameters may reveal a particular respiratory pathology, as exposed in [13]. A fractional order thermal model is also available in the literature [14]. This model arbitrarily proposes a Havriliak-Negami transfer function which relates the internal temperatures. Even though the fractional order model proves to be more accurate than an integer order one, the black box approach implies a lack of physical meaning for the model parameters.

Living tissue is normally perfused by blood, which is a natural thermal regulating mechanism. Perfusion allows tissue nutrition and gas exchange as well. A modified heat equation was proposed by H.H. Pennes in [15] in order to take into account additional phenomena:

$$ \rho c \frac{\partial T}{\partial t} = k \nabla^2 T + \rho_b c_b \omega_b (T_a - T) + q_m + q_{ext} \quad (1) $$

where $\rho_b c_b \omega_b (T_a - T)$ is the perfusion term, $q_m$ is metabolic heat generation and $q_{ext}$ is an external heat source.

Besides thermal modeling, this paper also deals with system identification. As an online system identification is sought, recursive methods are preferred. Recursive least squares with State Variable filter as well as Recursive Instrumental Variable methods have already been studied and applied for integer LPV systems [16]. For fractional order systems, one may find recursive methods able to recursively estimate coefficients and differentiation orders [17]. These methods require the calculation of an estimated output $\hat{y}(t)$, which is often achieved through an auxiliary model. As a consequence, calculation time may be constantly increasing as iterations progress. A simple truncated method of calculation may be used in order to have a fully functional identification method [18].

The paper is organized as follows: section II presents thermal two-port network formalism as well as its modifications in order to include blood perfusion. An analysis of the equivalent thermal impedance is presented in section III, system identification simulations are presented in section IV and final conclusions are given in section V.

## II. TWO-PORT NETWORK FORMALISM

Consider a thermal wall with 1D heat conduction as shown in figure 1, with section $S_w$ and length $L$. Its properties are given by conductivity $k$, density $\rho$ and specific heat capacity $c$. The blood thermal properties are given by $\rho_b$, $c_b$ and the perfusion coefficient is $\omega_b$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. 1D thermal system</figcaption></figure>

### A. Classic two-port network

For this section, physiological effects of blood perfusion are neglected. The system is analyzed as a continuous media with properties $k$, $\rho$ and $c$.

Starting with the heat equation:

$$ \rho c \frac{\partial T}{\partial t} = k \nabla^2 T, \quad (2) $$

and as conduction only goes in $x$ direction, it can be written as:

$$ \rho c \frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial x^2}. \quad (3) $$

The Laplace transform of equation (3), with null initial conditions, leads to:

$$ sT(x, s) = \frac{k \partial^2 T(x, s)}{\partial x^2}, \quad (4) $$

where $s$ denotes the Laplace variable.

By considering a heat flux input at $x = 0$ and a heat flux output at $x = L$, it then comes [19]:

$$ \phi_{in}(s) = -k S_w \frac{\partial T(x,s)}{\partial x} \Big|_{x=0} \quad (5) \\ \phi_{out}(s) = -k S_w \frac{\partial T(x,s)}{\partial x} \Big|_{x=L} $$
or even, put under a matrix form:

$$ \begin{bmatrix} T_{in}(s) \\ \phi_{in}(s) \end{bmatrix} = \mathbf{M} \begin{bmatrix} T_{out}(s) \\ \phi_{out}(s) \end{bmatrix} \qquad (6) $$

where

$$ \mathbf{M} = \begin{bmatrix} \cosh(\delta L) & \frac{1}{kS_w\delta}\sinh(\delta L) \\ kS_w\sinh(\delta L) & \cosh(\delta L) \end{bmatrix} \qquad (7) $$

with $\delta = \sqrt{\frac{s}{a}}$ and $a = \frac{\rho c}{k}$.

Note that $\delta$ involves the presence of a fractional (non-integer) operator, as its time-domain representation is a half-order derivative.

## B. Bio-heat two-port network

If the effect of blood perfusion is now to be included, one may take Pennes' equation (neglecting terms $q_m$ and $q_{ext}$):

$$ \rho c \frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial x^2} + \rho_b c_b \omega_b (T_a - T) \qquad (8) $$

Internal organ temperature is usually close to a constant core temperature $T_b$. Therefore, a variable change is introduced:

$$ \tilde{T} = T - T_b \qquad (9) $$

As $T_b$ is assumed constant, time and space derivatives for $\tilde{T}$ and $T$ will be equivalent. On the other hand, under a normal situation, blood tends to have a temperature close to the core:

$$ T_a \approx T_b \qquad (10) $$

This leads to the following:

$$ \rho c \frac{\partial \tilde{T}}{\partial t} = k \frac{\partial^2 \tilde{T}}{\partial x^2} - \rho_b c_b \omega_b \tilde{T} \qquad (11) $$

By defining thermal diffusivity $a$ and a perfusion heat coefficient $h_b$:

$$ a = \frac{k}{\rho c}, \quad h_b = \frac{\rho_b c_b \omega_b}{k} \qquad (12) $$

Equation (11) becomes:

$$ \frac{1}{a} \frac{\partial \tilde{T}}{\partial t} = \frac{\partial^2 \tilde{T}}{\partial x^2} - h_b \tilde{T} \qquad (13) $$

By taking the Laplace transform and assuming null initial conditions:

$$ \left[ \frac{s}{a} + h_b \right] \tilde{T}(x, s) = \frac{d^2 \tilde{T}(x, s)}{dx^2} \qquad (14) $$

This leads to an extremely similar two-port network structure as obtained in the previous subsection. There is however a difference for the $\delta$ definition:

$$ \delta_{bio} = \sqrt{\frac{s}{a} + h_b} \qquad (15) $$

## C. Circuit representation

Matrix $\mathbf{M}$ obtained for both the classical and bio-heat two-port network leads to a T circuit representation of the system, as shown in figure 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Thermal two-port network</figcaption></figure>
Impedance are given by:

$$Z_1(s) = Z_2(s) = \frac{1}{kS_w\delta}[\coth(\delta L) - \censch(\delta L)] \quad (16)$$

and

$$Z_3(s) = \frac{1}{kS_w\delta}\censch(\delta L) \quad (17)$$

where coth and csch respectively are the hyperbolic cotangent and cosecant functions. Variable $\delta$ will be taken as $\sqrt{\frac{s}{a}}$ for systems modelled after classic heat equation and as $\delta_{bio}$ if a model based upon bio-heat equation is used.

### D. Global Lung Model

A lung scheme is provided in figure 3. Each branch $n$ has a length $L_n$ as well as its cross-section $S_n$. A combination of Murray's law with a mechanical impedance model [20] allows to establish a relationship between consecutive cross-sections as well as lengths:

$$\frac{S_{n+1}}{S_n} = 2^{-\frac{2}{3}}, \quad \frac{L_{n+1}}{L_n} = 2^{-\frac{1}{3}}. \quad (18)$$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Lung scheme</figcaption></figure>

A symmetric lung model may be represented as the circuit shown in figure 4. On a same branching level $n$, the equivalent impedance is the parallel combination of $Z_1$ and $Z_3$, each lung level having $2^n$ branches which give:

$$Z_{i-eq,n}(s) = \frac{Z_{i,n}(s)}{2^n} \quad i = 1, 3. \quad n = 0, 1, \dots, 23. \quad (19)$$

### III. EQUIVALENT THERMAL IMPEDANCE

Lungs are usually perfused in a non-uniform way. Entrance airways are usually poorly perfused and act mostly as a passive tube system for air transportation. Gas exchange zone begins around $n = 14$. Therefore, a combination of the classic thermal two-port network and the bio-heat one is proposed to obtain a global lung model. This model is called *Mix* model as it mixes both two-port networks. The *Mix* model structure as well as the simulation parameters are shown in tables I and II, respectively. In table II, subscript *a* stands for air, whereas *t* stands for tissue.

<table><caption>TABLE I
*Mix* MODEL STRUCTURE</caption><thead><tr><th>Model</th><th>Dead space (n &lt; 14)</th><th>Exchange space (n ≥ 14)</th></tr></thead><tbody><tr><td>Mix</td><td>Classic</td><td>Bio-heat</td></tr></tbody></table>

Let us consider the entrance level for the perfused zone $n = 14$. The thermal equivalent impedance results from the complex expressions in (16) and (17), as well as the complex circuit network shown in figure 4. In order to approximate the exact equivalent thermal impedance of the lungs, a simple first kind fractional order model is proposed:

$$Z_{\text{lung},n}(s) = \frac{K}{1 + \left(\frac{s}{\omega_c}\right)^{\nu}}. \quad (20)$$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Lung equivalent circuit model</figcaption></figure>

<table><caption>TABLE II
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>

The general form of this model is:

$$G(s) = \frac{K}{\tau s^{\nu} + 1} \qquad (21)$$

By taking the static gain from the exact frequency response and performing a simple optimization with FPA algorithm for parameter $\tau$ [21], one obtains the following approximated model:

$$G_{app-lung}(s) = \frac{3.48}{8.794s^{0.5} + 1} \qquad (22)$$

Gain diagrams for both the exact and the approximated frequency responses are given in figure 5.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Exact and approximated input equivalent impedance at n = 14</figcaption></figure>

## IV. SYSTEM IDENTIFICATION

For system identification simulations, a recursive two-stage algorithm will be used in order to estimate both the coefficients and the differentiation order. The used method is the LMRPEM-2 which was successfully validated in [17].
Classic prediction error method minimizes a quadratic error criterion:

$$J = \frac{1}{2} \sum_{k=1}^{t} \epsilon(k)^2 \qquad (23)$$

where the error function is originally defined as:

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

where $y^*(t)$ stands for the noisy output and the estimated output $\hat{y}(t)$ is computed as:

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

However, the nature of fractional order derivatives lead to the choice of an extended error in order to take into account fractional system natural long-memory property.

The error signal $\tilde{\epsilon}$

$$\tilde{\epsilon}(k) = [\epsilon(0), \epsilon(T_s), \epsilon(2T_s), \dots, \epsilon(kT_s)]^T \qquad (26)$$

gathers the errors at all instants from $t = 0$ to the current time $t = kT_s$ and is increased by one data-point per iteration.

In this method, the extended error is used and at each single iteration, updates coefficients as well as updated differentiation orders are estimated. For further information, the complete description of the algorithm may be available in [17].

For simulations, the sampling time is set to $T_s = 1s$. A multisinus heat density flux signal is used as an input with an amplitude going from $-10$ to $10mW.m^{-2}$. This signal contains 10 harmonics. The output is calculated by using the lung approximated model shown in previous sub-section. An $SNR = 20dB$ is used for this simulated data. The estimated coefficients are initialized at $\hat{K}_0 = 30$, $\hat{\tau}_0 = 10$ and $\hat{\nu}_0 = 0.5$. The coefficients for LMRPEM-2 algorithm are set to $\gamma_\rho = 0.05$ and $\gamma_\nu = 0.05$. Input-output data are shown in figure 6.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Input-output data with $SNR = 20dB$</figcaption></figure>

### A. Parameter estimation without truncation

For this subsection, LMRPEM-2 method will be applied as described in [17]. Figure 7 shows the result for both coefficients and differentiation order estimation.

As can be observed, there is an important transient at the start of the estimation regarding $\tau$ and $\nu$. This may be produced by an insufficient information in order to correctly determine the system dynamics (which is slow in this example). After this transient phase, there is a convergence towards the true values of the parameters. By the end of the iterations, the estimated model is given by:

$$\hat{G}_{lung}(s) = \frac{3.47}{8.80s^{0.5} + 1}, \qquad (27)$$

which is very close to the true model (22).

Even though this estimation looks promising, there is a main implementation issue when applying LMRPEM-2 method. Each iteration requires the calculation of the estimated output $\hat{y}(t)$, which is achieved through a convolution product between the input and an auxiliary model. As iterations progress, the calculation of this output becomes progressively slower. This, of course, implies that the calculation time increases as iterations do. Figure 8 shows the calculation time for the previous identification.
<figure><img src="image_8.png" /><figcaption>Fig. 7. LMRPEM-2 estimation without truncation</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Calculation time $t_{calc}$ vs iterations with the LMRPEM-2 without truncation</figcaption></figure>

As one observes, there is a relatively fast increase in calculation time as iterations progress. By the end of the iterations, the calculation time for a single iteration is around 1s, which is the sampling time. If the calculation time increases beyond the value of the sampling time, the algorithm will no longer be exploitable in real-time scenarios in spite of its recursive formulation. As a consequence, a simple truncated method may be used in order to limit the constant increase in the calculation time.

## B. Parameter estimation with truncation

In order to limit the calculation time growth, a truncated approximated method is proposed: the convolution to estimate $\hat{y}$ will be taken over a finite temporal length $L$ instead of using the whole past of the input. At a given moment $t^{act}$, the output will be approximated by a convolution between $t = t^{act} - L$ and $t = t^{act}$. Even though this introduces some error due to inaccurate initial conditions, it may lead to adequate parameter estimation for stable systems and a long enough window $L$. The whole developments on how to fix a suitable memory length $L$ to limit the approximation errors are proposed in [18].

In this example, the truncated LMRPEM-2 algorithm is applied by applying a memory window length of $L = 1500$s. Figure 9 shows parameter estimation evolution and the calculation time is shown in figure 10.

The parameter estimation is still accurate, even though the use of a truncation leads to small oscillations for parameter $τ$ (barely visible). On the other hand, the calculation time at each iteration has effectively been limited and is under the sampling time of $T_s = 1$s for all the iterations. Therefore, the algorithm is now fully functional for a real-time application.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Truncated LMRPEM-2 estimation with truncation length $L = 1500s$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Calculation time $t_{calc}$ vs iterations with the truncated LMRPEM-2</figcaption></figure>

V. CONCLUSIONS

A two-port network formalism is applied in order to obtain accurate thermal circuit models. This model is extended by combining it with the bio-heat equation in order to take into account the blood perfusion role on thermal regulation, which may be important in physiological scenarios. A global circuit model combining classic and bio-heat circuits is proposed and its frequency response is approximated by a simple first kind fractional order model. The second part of the paper considers system identification using the simple first kind model obtained from the global lung thermal circuit. An innovative recursive two-stage algorithm suited for fractional order identification is applied and both coefficients as well as differentiation orders are successfully estimated. Finally, the importance of limiting the calculation time at each iteration is exposed and a simple truncated method, the LMRPEM-2 one, is applied in order to guarantee real-time system identification with fractional models. Further studies could consider cylindrical and spherical geometries for the lung thermal model and other equivalent thermal impedance expressions may also be considered.

REFERENCES

[1] R. Badenes, A. Lozano, and F. J. Belda, "Postoperative Pulmonary Dysfunction and Mechanical Ventilation in Cardiac Surgery," Feb. 2015. [Online]. Available: https://www.hindawi.com/journals/ccrp/2015/420513/

[2] R. E. Giesinger and P. J. McNamara, "The Impact of Therapeutic Hypothermia on Pulmonary Hemodynamics of Meconium Aspiration Syndrome," *American Journal of Respiratory and Critical Care Medicine*, vol. 198, no. 2, pp. 286-287, Mar. 2018, publisher: American Thoracic Society - AJRCCM.

[3] J. Hildebrandt, "Dynamic properties of air-filled excised cat lung determined by liquid plethysmograph." *Journal of Applied Physiology*, vol. 27, no. 2, pp. 246-250, 1969.
[4] —, "Pressure-volume data of cat lung interpreted by a plastoelastic, linear viscoelastic model," *Journal of Applied Physiology*, vol. 28, no. 3, pp. 365-372, 1970.

[5] A. B. DuBois, A. W. Brody, D. H. Lewis, and B. F. Burgess, "Oscillation mechanics of lungs and chest in man," *Journal of Applied Physiology*, vol. 8, no. 6, pp. 587-594, 1956.

[6] D. Navajas, R. Farre, J. Canet, M. Rotger, and J. Sanchis, "Respiratory input impedance in anesthetized paralyzed patients," *Journal of Applied Physiology*, vol. 69, no. 4, pp. 1372-1379, 1990.

[7] B. Diong, H. Nazeran, P. A. Nava, and M. Goldman, "Modeling human respiratory impedance," *IEEE Engineering in Medicine and Biology Magazine*, vol. 26, pp. 48-55, 2007.

[8] E. R. Weibel, "Mandelbrot's fractals and the geometry of life: A tribute to benoît mandelbrot on his 80th birthday," *Fourth International Symposium on Fractals in Biology and Medicine*, 2014.

[9] C. Ionescu, R. De Keyser, J. Sabatier, A. Oustaloup, and F. Levron, "Low frequency constant-phase behavior in the respiratory impedance," *Biomedical Signal Processing and Control*, vol. 6, no. 2, pp. 197-208, 2011.

[10] C. Ionescu, "Relationship between anatomical structure, fractional order models and fractal dimension in the respiratory impedance of healthy patients," *07 2012*, pp. 1425-1430.

[11] M. Nakagawa and K. Sorimachi, "Basic characteristic of fractance device," *IEICE Trans. Fund*, 1992.

[12] B. Krishna and K. Reddy, "Active and passive realization of fractance device of order 1/2," *Active Passive Electronic components*, 2008.

[13] D. Copot, "Fractional calculus based methods and models to characterize diffusion in the human body," Ph.D. dissertation, University of Gent, 2018.

[14] M. Pellet, P. Melchior, Y. Abdelmoumen, and A. Oustaloup, "Fractional thermal model of the lungs using havriliak-negami function," vol. 3, 2011.

[15] H. H. Pennes, "Analysis of tissue and arterial blood temperatures in the resting human forearm," *Journal of Applied Physiology*, vol. 1, no. 2, 1948. [Online]. Available: https://journals.physiology.org/doi/abs/10.1152/jappl.1948.1.2.93

[16] A. Padilla, "Recursive identification of continuous-time systems with time-varying parameters," Ph.D. dissertation, Université de Lorraine, 2017.

[17] S. Victor, J.-F. Duhé, P. Melchior, Y. Abdelmounen, and F. Roubertie, "Long-memory recursive prediction error method for identification of continuous-time fractional models," *Nonlinear Dynamics*, vol. 110, no. 3, pp. 635-648, 2022.

[18] J.-F. Duhé, S. Victor, P. Melchior, Y. Abdelmounen, and F. Roubertie, "Fractional derivative truncation approximation for real-time applications," *Communications in Nonlinear Science and Numerical Simulation*, vol. 119, p. 107096, 2023. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S100757042300014X

[19] J. Duhé, 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.

[20] F. Kuwahara, Y. Sano, J. Liu, and A. Nakayama, "A porous media approach for bifurcating flow and mass transfer in a human lung," *Journal of Heat Transfer*, vol. 131, no. 101013, Jul. 2009. [Online]. Available: https://doi.org/10.1115/1.3180699

[21] X. Yang, "Flower pollination algorithm for global optimization," *Unconventional Computation and Natural Computation*, 2012.