# Lung Thermal Transfer System Identification With Fractional Models

Stephane Victor, Pierre Melchior, Mathieu Pellet, Alain Oustaloup

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

## To cite this version:

Stephane Victor, Pierre Melchior, Mathieu Pellet, Alain Oustaloup. Lung Thermal Transfer System Identification With Fractional Models. IEEE Transactions on Control Systems Technology, 2020, 28 (1), pp.1-11. 10.1109/TCST.2018.2877606. hal-02010639

HAL Id: hal-02010639

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

Submitted on 9 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Lung thermal transfer system identification with fractional models

Stéphane Victor, Pierre Melchior, Mathieu Pellet and Alain Oustaloup

## Abstract

System identification through fractional models is proposed for modeling thermal transfers in the lungs. In general, during open-heart surgery, an extracorporeal circulation (ECC) is carried out on the patient. First, the lungs are disconnected from the circulatory system so that an artificial heart/lung machine can be plugged on the blood stream. Pulmonary cells die when they are no more supplied with fresh blood which may result in postoperative respiratory complications for the patient. Also, bronchial hypothermia is carried out to protect the lungs. Indeed, cooling the lung organ helps slowing down its deterioration. Unfortunately the lung thermal properties are not well-known yet, therefore mathematical models are useful and needed in order to improve the knowledge of these organs. As proved by several previous works, fractional models are well-suited to model the dynamics of fractal systems and to model thermal systems thanks to its model parameter compactness but also thanks to the solving of the heat equation whose analytical result reveals a fractional differentiation operator of order 0.5. System identification results are provided on real data measured on sheep lungs. Thus, this paper studies the comparison between a rational model and two kinds of fractional models, a classic one and another one using the Havriliak-Negami function. Time domain and frequency domain analysis are provided showing the efficiency of fractional models, but more precisely of Havriliak-Negami transfer functions.

## Index Terms

System identification, Bio-engineering, Fractional order model, lung modeling, non-integer differentiation, black-box identification

## I. INTRODUCTION

NOWING the thermal behavior of lungs is crucial for anesthesiologists and heart cardiologist surgeons. Also, an original study for thermal modeling of lungs while open-heart surgery is proposed. In most open-heart surgery operations, it is necessary to stop the heart from beating so that the surgeon may operate properly. An extracorporeal circulation (ECC) is then provided thanks to an artificial heart/lung device connected to the circulatory system of the patient. Such a device does ensure the organ functions, namely heart and lung functions, however, the organs themself are therefore excluded from the blood stream and consequently deprived of blood and oxygen during the whole surgery, leading to inevitable damages for the organ tissues recognized as postoperative pulmonary dysfunctions. These damages may be reversible for younger patients but can cause major complications for elder ones. Anesthetists have developed methods to protect the lungs and to restrain the damages through hypothermia, that is to say by cooling them. Hypothermia enables slowing the metabolism of the lung cells and therefore slowing damage progression, for example, by directly applying ice bags or fresh gas intern insufflation. These techniques are unfortunately empiric, also, there is a need for real physical models that well describes the lung temperature evolution during hypothermia. Moreover, fresh gas insufflation has good results, nevertheless, it is based on the assumption that the lung thermal insulation increases during ECC, an assumption that has not been verified yet. A lung thermal model during ECC and during standard conditions will help to verify this assumption. Different models of thermal transfers are identified and presented in this paper from real data measured from sheep lungs.

In biology, lungs can be represented like fractal geometry systems. As it is known, non-integer differentiation is closely linked to recursive and fractal geometry: in mechanics such as a porous dyke [1], thermal diffusive phenomena [2][3], charge diffusion in acid batteries [4], in rheology [5], Foucault currents inside rotor bars in induction machines [6], in natural environment such as coastal geometry [7], in biology systems [8, 9, 10, 11]... Left and right lungs being not exactly symmetrical, their geometry is nevertheless not strictly fractal. However each lung, left or right, does have a fractal geometry [9][12, chap. 6]. Therefore, modeling a lung with fractional differentiation is completely justified [9, 13, 14][12, chap. 6]. Research works using fractional differentiation as modeling tool of the lungs has been led notably to model the lung respiratory function and to study the link between the impedance variations and eventually respiratory illnesses [15]. In this paper, the works hereby presented have been carried out in the Anesthesiology - Intensive care service of the Haut-Lévêque Hospital of Bordeaux, and more specifically, with Dr Youssef Abdelmoumen. This service, specialized in open-heart surgery has developed a fresh air insufflating technique for bronchial hypothermia. The improvement of this technique requires mathematical models with parsimonious parameters, easily exploitable and rending the thermal changes inside of the lung during surgery that needs extracorporeal circulation.

Stéphane Victor, Pierre Melchior, Mathieu Pellet, and Alain Oustaloup are researchers in IMS, UMR 5218 CNRS - Université de Bordeaux - Bordeaux INP 351 cours de la Libération, 33405 TALENCE cedex, FRANCE (e-mail: firstname.lastname@ims-bordeaux.fr)

This project has been carried out with the collaboration of the IMS laboratory and the Anesthesiology - Intensive care service of the Haut-Lévêque Hospital of Bordeaux.

We would like to thank Dr. Youssef Abdelmoumen, surgeon at the Anesthesiology - Intensive care service of the Haut-Lévêque Hospital of Bordeaux. Manuscript received February 28, 2018.
In order to obtain such models, several measures have been achieved by the medical team of Haut-Lévêque hospital on sheep. After presenting the experimental protocol, measures have been used for system identification with a black box approach. In this context of lung fractality, thermal system and parameter parsimony, non-integer models are particularly suited, and especially the Havriliak-Negami function that has shown promising results in system identification in terms of parameter parsimony [16]. Preliminary studies have been carried out with lung thermal system identification with fractional models by using a simple first order fractional model called as explicit form<fnref n="1" /> models [13]: $M(s) = \frac{a_0 s^{\nu_1}}{1 + b_1 s^{\nu_2}}$. Moreover, another preliminary study has compared Havriliak-Negami models with a simple rational transfer function of the form [17]: $M^{rat}(s) = \frac{\sum_{i=0}^{n} b_i}{1 + \sum_{j=1}^{4} f_j s^j}$. The promising results have encouraged the authors to continue their research work in system identification with Havriliak-Negami transfer functions. Thus, the last section of this paper presents the latest improved results with a comparison between a non-integer model of Havriliak-Negami type and a more complete rational one on the one hand, and between a non-integer one and a more complete fractional approximation on the other hand. Furthermore, a frequency domain analysis is also provided at the end of this paper. All these results can be found in [18].

Section II relates the modeling of the lung respiratory function from its mechanical properties. Then after presenting the experimental protocol, system identification is proposed for temperature transfers in section III. Section IV presents the experimental results with several model performance comparisons both in time-domain and frequency domain to finally conclude in section V.

## II. LUNG RESPIRATORY FUNCTION MODELING

Authors in [15] have tried to diagnosis lung respiratory illnesses with fractional models. Measures were confronted with the fractional model proposed by [12, chap 6]: the lung input impedance is defined as the ratio between the pressure and the air flow entering the lung. The idea of this modeling relies on the electrical analogy of the mechanical properties of the lung. Indeed, the respiratory duct can be represented by an arborescent network of 24 levels. In the respiratory viewpoint, a lung can be modeled by a bifurcation arborescent structure where each branch presents (in electrical equivalence) (see Figure 1):

* a dissipative effect which can be characterized by a resistance ($R_i$), due to viscosity and turbulence energy losses;

* an inductive effect which can be characterized by an inductance ($L_i$), due to inertial effects of motion air;

* a capacitive effect which can be characterized by a capacitance ($C_i$), due to the conjoint elasticity of the air and the branch which makes (by a tank effect) the ingoing flow different from the outgoing one.

<figure><img src="image_3.png" /><figcaption>Fig. 1: Electrical analogy: arborescent arrangement of RLC cells</figcaption></figure>

On Fig. 1, The current $I_0$ and the voltage $U_0$ are the images of, respectively, the pressure and the air flow entering the lung. Moreover, the impedances are expressed by:

$$z_i(j\omega) = R_i + jL_i\omega.$$

By gathering in one single node the node set of same potential, and by repeating this operation for each network potential, the arborescent arrangement of gamma RLC cells is brought to a cascade arrangement of gamma RLC cells (Fig. 2) in which:

$$Z_i(j\omega) = \frac{z_i(j\omega)}{2^{i-1}} = R'_i + jL'_i j\omega,$$

with $R'_i = R_i/2^{i-1}$, $L'_i = L_i/2^{i-1}$ and $C'_i = 2^{i-1}C_i$.

<footnote n="1">¹A fractional model of the explicit form has the following form: $M^{frac}(s) = \frac{\sum_{i=0}^{n} b_i s^{\nu_i}}{\sum_{j=0}^{m} a_j s^{\nu_j}}$ where the coefficients $a_i, b_j \in \mathbb{R}$ and the fractional differentiation orders $\nu_i, \nu_j \in \mathbb{R}^+$, s being the Laplace variable</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 2: Cascade arrangement of gamma RLC cells</figcaption></figure>

Under the assumption of the resistance, inductance and capacitance recursivity of the arborescent network of Fig. 1, linked to the fractal geometry of the lung, namely:

$$R_{i+1} = \alpha R_i, L_{i+1} = \lambda L_i \text{ and } C_{i+1} = \eta C_i,$$

or even,

$$R'_{i+1} = \alpha' R'_i, L'_{i+1} = \lambda' L'_i \text{ and } C'_{i+1} = \eta' C'_i,$$

with $\alpha' = \alpha/2$, $\lambda' = \lambda/2$ and $\eta' = 2\eta$, the total admittance can be expressed by ($s$ being the Laplace variable):

$$Y(s) = \frac{\frac{1}{R'_{1} + L'_{1}s}}{1 + \frac{\frac{1}{[C'_{1}s(R'_{1} + L'_{1}s)]}}{\frac{1}{[C'_{2}s(R'_{2} + L'_{2}s)]}} \cdot \frac{1 + \dots + \frac{1}{[C'_{N}s(R'_{N} + L'_{N}s)]}}{1 + \frac{1}{1 + [C'_{N}s(R'_{N} + L'_{N}s)]}}}$$

After some changes of variable (see [12, chap. 6]), this admittance can be rewritten into:

$$Y(s) = \frac{1/R_1}{K(\alpha, \eta)} (R_1 C_1 s)^n,$$

a form that reveals a fractional order differentiation such that

$$n = \frac{\log \alpha}{\log(\alpha\eta)}$$

and $K(\alpha, \eta)$ is defined such that $K(\alpha, \eta)K(\eta, \alpha) = \alpha^n$.

The analytical expression of the arborescent network impedance can thus be defined as a fractional transfer function where a physical interpretation can be given to the parameters. Moreover, only five parameters need to be identified, namely, $R_1$, $C_1$ and the recursive factors. This type of model also presents the advantage to rely on an electrical analogy, and has therefore a generic character: a similar modeling process can be applied to model physical phenomena other than the respiratory system as long as it is possible to define an effort/flux relation in the Bond Graph sense, analogous to the voltage/current relation proposed here.

### III. SYSTEM IDENTIFICATION

<figure><img src="image_5.png" /><figcaption>Fig. 3: Experimental protocol: anesthetized sheep equipped with sensors (left) and artificial heart/lung device with ECC start (right)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4: Rough diagram of lungs (• denoting the approximate location of the thermal sensors)</figcaption></figure>

A. Experimental protocol

Dr Youssef Abdelmoumen and his team are specialized in open-heart surgery. They provided the thermal data acquired on sheep. The experiment was carried out at the PTIB platform of Haut-Lévêque Hospital.

Lungs are organs presenting an arborescent (or tree-like) structure. From the trachea, the respiratory duct is subdivided in series of bronchial tubes that becomes more and more narrow. The general aim of this work is to study the variation of thermal characteristics of the lungs during extracorporeal circulation on a sheep. Thus, the medical team has recorded temperatures of the lungs during open-heart surgery before the ECC and during the ECC. The medical team has then set up an experimental protocol enabling to measure the temperature at different levels of the lung tree structure inside bronchus. Animal testing are led on sheep, the animals being taken care by experimented surgeons and in conformity with the Animal Care and Use Committee. The sheep is first anesthetized, then intubated and the different sensors necessary to the measures are introduced inside its lungs via the trachea. An opening is then practiced at the sternum level enabling to achieve the different operations necessary for the ECC connection and start (see Fig. 3).

All temperatures are measured thanks to GI type thermocouple. These thermocouples are linked to a National Instruments PXI interface. A special software was designed with LABVIEW® and load into the PXI in order to acquire and record the signals at a sampling frequency of 50Hz.

Temperatures are recorded thanks to three sensors placed at the entry of the three first bronchus and a fourth sensor is used to measure the body temperature (see Fig. 4). Measures are limited to these three locations because of the difficulty to accurately plant a sensor deeper into the lung. During the experiment, cold and hot flow are alternated at the trachea intubation.

B. Temperature transfers

The temperatures measured on the bronchus $T_1$, $T_2$ and $T_3$ are normalized by the sheep body temperature $T_c$. From the measured data, two temperature transfers can be defined: one from bronchus 1 to bronchus 2

$$H_{12}(s) = \frac{\Delta T_2(s)}{\Delta T_1(s)} \qquad (1)$$

and one from bronchus 2 to bronchus 3

$$H_{23}(s) = \frac{\Delta T_3(s)}{\Delta T_2(s)}, \qquad (2)$$

where

$$\Delta T_k(s) = T_k(s) - T_c(s) \quad \text{for } k = 1, 2, 3.$$

Putting into equation these transfers turns out to be tricky. Indeed, the multiplicity of conduction/convection phenomena and the ignorance of the thermal properties of the lung biologic tissues prevent to determine an analytical model. Therefore, a black-box model is proposed for system identification.
Non-integer differentiation has enabled to extend the possibilities of time and frequency domain system identification beyond the integer case for differentiation order. Cole-Cole and Davidson-Cole transfer functions helps to model complex diffusive systems such as thermal systems with few parameters (parameter parsimony) [19]. For similar modeling error, such systems represented by rational transfer functions of integer order would need more than twice the number of parameters used for fractional transfer functions [20].

Furthermore, the Havriliak-Negami function provides another freedom degree in system identification by simultaneously introducing a fractional character in intern and extern orders of Cole-Cole and Davidson-Cole transfer functions. Also, the Havriliak-Negami function is particularly suited in this black-box system identification in order to model the transfers $H_{12}(s)$ and $H_{23}(s)$.

The only minor difficulty in the use of the Havriliak-Negami function lies on its approximation and its time simulation. The method used is the one proposed in [16] which is recalled in the following paragraph.

### C. The Havriliak-Negami function

Historically, the Havriliak-Negami function has been introduced to understand the frequency behavior of dielectric relaxation processes in some polymers [21]. Its definition

$$F_{HN}(s) = \frac{1}{\left(\left(\frac{s}{\lambda}\right)^{\nu_1} + 1\right)^{\nu_2}} \quad (3)$$

includes the Cole-Cole function definition, such as:

$$F_{CC}(s) = \frac{1}{\left(\frac{s}{\lambda}\right)^{\nu_1} + 1},$$

and the Davidson-Cole function one, such as

$$F_{DC}(s) = \frac{1}{\left(\frac{s}{\lambda} + 1\right)^{\nu_2}}.$$

The idea of this approximation rests on a decomposition of the frequency behavior of the Havriliak-Negami function into two parts: its asymptotic frequency behavior and the difference between its real and asymptotic behavior. In the same way as classic non-integer transfer functions can be approximated by a recursive distribution of poles and zeros, the difference between the real and asymptotic behavior can be synthesized with a particular disposition of poles and zeros.

First, the methodology is presented for $\nu_1 \in ]0, 1[$ as $s^{\nu_1} + \lambda$ has no root, then it will be generalized to $\nu_1 \in \mathbb{R}$ by taking into account of roots when $\nu_1 > 1$.

1) Synthesis for $0 < \nu_1 < 1$: Consider the Havriliak-Negami function $F_{HN}(s)$ given by:

$$F_{HN}(s) = (s^{\nu_1} + \lambda)^{\nu_2} \quad (4)$$

with $\nu_1 \in ]0, 1[$, $\nu_2 \in \mathbb{R}$ and $\lambda \in \mathbb{R}^+$. This function can be written under the following form:

$$F_{HN}(s) = F_1(s)\kappa(s) \quad (5)$$

in which

$$F_1(s) = (s + \lambda^{\nu_1})^{\nu_1 \nu_2} \quad (6)$$

is a Davidson-Cole function that enables considering the gain and phase asymptotic diagrams of the Havriliak-Negami function $F_{HN}(j\omega)$. Taking into account these asymptotic behaviors thus brings in a recursive distribution of poles and zeros (Oustaloup approximation of Davidson-Cole function).

Consequently, $\kappa(s)$ describes the term related to the difference between the real and asymptotic frequency behavior and is defined by:

$$\kappa(s) = \frac{F_{HN}(s)}{F_1(s)}, \quad (7)$$

and can be approximated by a non recursive distribution of poles and zeros.

Indeed, the gain diagram of function $\kappa(j\omega)$ shows a symmetry with respect to the transitional frequency $\lambda$. The synthesis of $\kappa(j\omega)$ then relies on as many poles and zeros distributions as frequency bands considered on either side of $\lambda$. This type of synthesis necessitates an identical number of frequency band for $\omega > \lambda$ and for $\omega \le \lambda$. This synthesized function $\tilde{\kappa}(j\omega)$ is thus based on a non recursive distribution of poles and zeros whose principle is presented on Figure 5.

The function $\kappa(j\omega)$ is synthesized on the frequency half-band $\Delta$ which is given by the following relations:

$$\begin{aligned} \Delta_1 &= \max\left(\frac{\lambda}{\chi \cdot \omega_h}, \frac{\chi \cdot \omega_h}{\lambda}\right) \\ \Delta_2 &= \max\left(\frac{\chi \cdot \lambda}{\omega_b}, \frac{\omega_b}{\chi \cdot \lambda}\right) \\ \Delta &= \max(\Delta_1, \Delta_2), \end{aligned} \quad (8)$$
<figure><img src="image_7.png" /><figcaption>Fig. 5: Principle of the non recursive pole (x) and zero (o) distribution for the synthesis of $\kappa(j\omega)$ on the frequency band $[\omega_b, \omega_h]$ when a) $\lambda \le \omega_b \le \omega_h$, b) $\lambda \le \omega_b \le \omega_h$, and c) $\omega_b \le \omega_h \le \lambda$</figcaption></figure>

where $\chi$ denotes the growth rate of the frequency band to take into account the edge effect ($\chi = 100$ limits significantly this effect).

The synthesized function $\tilde{\kappa}(s)$

$$ \tilde{\kappa}(s) = \prod_{k=1}^{K} \prod_{m=1}^{M_k} \frac{s + \omega_{k,m}'}{s + \omega_{k,m}} $$

relies on the poles and zeros of each of the K frequency bands $[\omega_{e_k}, \omega_{e_{k+1}}]$ such that

$$ \omega_{e_k} = \omega_b \delta^{k-1} \quad \text{with} \quad k = 1, \dots, K+1 $$

whose distribution is given by the relations identical to the ones used for the synthesis of a fractional differentiator, by
considering poles and zeros on frequency bands:

$$ \omega'_{k,0} = (\alpha_k)^{\frac{1}{2}} \eta_k \omega_{e_k} $$

$$ \omega_{k,0} = (\alpha_k)^{\frac{1}{2}} \omega_{e_k} $$

$$ \frac{\omega'_{k,m+1}}{\omega'_{k,m}} = \frac{\omega_{k,m+1}}{\omega_{k,m}} = \alpha_k \eta_k > 1 $$

$$ \frac{\omega_{k,m+1}}{\omega'_{k,m}} = \alpha_k > 0 $$

$$ \frac{\omega'_{k,m+1}}{\omega_{k,m}} = \eta_k > 0 $$

$$ \nu_k = \frac{\log \alpha_k}{\log(\alpha_k \eta_k)}. $$

The succession of poles and zeros thus established enables writing $\tilde{\kappa}(s)$ under the form:

$$ \tilde{\kappa}(s) = \prod_{k=1}^{K} \prod_{m=1}^{M_k} \frac{s + \omega'_{k,m}}{s + \omega_{k,m}} = \prod_{i=1}^{K \cdot M_k} \frac{s + \omega'_{i}}{s + \omega_{i}} $$

2) Synthesis for $\nu_1 \ge 1$: The value of order $\nu_1$ determines the number of $s$ poles of $s^{\nu_1} + \lambda$. The definition of the Havriliak-Negami function $(s^{\nu_1} + \lambda)^{\nu_2}$ introduces a $(s - s_m)^{\nu_2}$ for each $s_m$ pole. Thus, the approximation of the Havriliak-Negami function involves $M$ eigenvalues such as:

$$ F_2(s) = \prod_{m=1}^{M} (s - s_m)^{\nu_2} $$

enabling writing the Havriliak-Negami function under the form:

$$ F_{HN}(s) = F_2(s) \cdot F_1(s) \cdot \kappa_2(s), $$

where $F_1(s)$ enables respecting the gain and phase asymptotic behaviors of the Havriliak-Negami function $F_{HN}(s)$:

$$ F_1(s) = (s + \lambda)^{\nu_2(\nu_1 - M)}. $$

The synthesis of the Havriliak-Negami function can then be put under the form:

$$ F_R(s) = \prod_{m=1}^{M} (s - s_m)^{\nu_2} \cdot (s + \lambda)^{\nu_2(\nu_1 - M)} \cdot \kappa_2(s). $$

Involving the approximations of the pole and zero distribution relative to each terms, the Havriliak-Negami function is then given by the relation:

$$ F_R(s) = \prod_{m=1}^{M} C_m \prod_{k_1=1}^{K_1} \frac{s + \omega'_{mk_1}}{s + \omega_{mk_1}} \cdot \prod_{k_2=1}^{K_2} C_2 \frac{s + \omega'_{k_2}}{s + \omega_{k_2}} \cdot \prod_{k_3=1}^{K_3} C_3 \frac{s + \omega'_{k_3}}{s + \omega_{k_3}} $$

where the frequencies $\{\omega_{mk_1}, \omega'_{mk_1}\}$ and $\{\omega_{k_2}, \omega'_{k_2}\}$ are given by similar relations to those leading to the synthesis of a non-integer differentiator ($k = mk_1$ or $k_2$):

$$ \omega_H = \omega_h \cdot \chi $$

$$ \omega'_0 = \alpha^{\frac{1}{2}} \eta \lambda $$

$$ \omega_0 = \alpha^{\frac{1}{2}} \lambda $$

$$ \frac{\omega'_{k+1}}{\omega'_{k}} = \frac{\omega_{k+1}}{\omega_{k}} = \alpha \eta > 1 $$

$$ \frac{\omega_{k+1}}{\omega'_{k}} = \alpha > 0 $$

$$ \frac{\omega'_{k+1}}{\omega_{k}} = \eta > 0. $$

(9)
$\kappa_2(s)$ is achieved on the frequency band $[\omega_{b_2}, \omega_{h_2}]$ by satisfying relations (8) and such that:

$$\omega_{h_2} = \lambda\Delta$$

$$\omega_{b_2} = \frac{\lambda}{\Delta}$$

$$\delta = \Delta^{\frac{1}{K}}.$$

The poles and zeros $\{\omega_{k_3}, \omega'_{k_3}\}$ of each of the K frequency bands $[\omega_{e_k}, \omega_{e_{k+1}}]$ such that:

$$\omega_{e_k} = \omega_B \delta^{k-1} \quad \text{with} \quad k = 1, \dots, K+1$$

are obtained by relations identical to those used for the synthesis of a non-integer differentiator by considering $M_k$ poles and zeros per frequency band (see relations in (9)).

3) Synthesis for $\nu_1 \le 0$: Now, the relation

$$F_{HN}(s) = (s^{\nu_1} + \lambda)^{\nu_2} \quad (10)$$

can be put under the form:

$$F_{HN}(s) = (s^{-|\nu_1|} + \lambda)^{\nu_2} = \lambda^{\nu_2} \left( \frac{1}{\lambda s^{|\nu_1|}} + 1 \right)^{\nu_2} \quad (11)$$

or even

$$F_{HN}(s) = \frac{\lambda^{\nu_2}}{s^{|\nu_1|\nu_2}} \left( \frac{1}{\lambda} + s^{|\nu_1|} \right)^{\nu_2}, \quad (12)$$

where the first term represents a fractional differentiator of order $|\nu_1|\nu_2$, and the second one being a Havriliak-Negami function dealt previously. The synthesis $F_R(s)$ of the function $F_{HN}(s)$ is then given by:

$$F_R(s) = \prod_{k_0=1}^{K_0} C_0 \frac{s + \omega'_{k_0}}{s + \omega_{k_0}} \cdot \prod_{m=1}^{M} C_m \prod_{k_1=1}^{K_1} \frac{s + \omega'_{mk_1}}{s + \omega_{mk_1}} \\ \cdot \prod_{k_2=1}^{K_2} C_2 \frac{s + \omega'_{k_2}}{s + \omega_{k_2}} \cdot \prod_{k_3=1}^{K_3} C_3 \frac{s + \omega'_{k_3}}{s + \omega_{k_3}}$$

where the poles and zeros $\{\omega_{k_0}, \omega'_{k_0}\}$ are defined by relations enabling the synthesis of a non-integer integrator.

This synthesis method enables approximating the Havriliak-Negami function on a frequency band given by the real orders $\nu_1$ and $\nu_2$ and thus simulating it easily for system identification purpose.

## IV. EXPERIMENTAL RESULTS

Fractional models are well-suited for modeling diffusive phenomena and therefore for thermal system modeling. In time domain system identification, the Havriliak-Negami function, still not much used, has shown interesting results for modeling thermal diffusion on a metallic rod [22] or even for modeling thermal transfers on high speed processing tool [16].

The Havriliak-Negami is first compared with a rational (or integer) transfer function to highlight the interest of non-integer models in this study context. Then, in order to measure the interest of the specific Havriliak-Negami function in system identification, the obtained model is compared with a fractional model.

Two kinds of data are identified: one during the normal operating of the organ before starting the ECC, and the other one during the ECC.

The temperatures measured in the bronchus $T_1, T_2$ and $T_3$ are normalized with the animal body temperature $T_c$. Two transfers are then defined: the transfer from to the first lung bronchus relative temperature $\Delta T_1$ towards the second bronchus relative temperature $\Delta T_2$

$$H_{12}(s) = \frac{\Delta T_2(s)}{\Delta T_1(s)} \quad (13)$$

and the transfer from to the second lung bronchus relative temperature $\Delta T_2$ towards the third bronchus relative temperature $\Delta T_3$

$$H_{23}(s) = \frac{\Delta T_3(s)}{\Delta T_2(s)} \quad (14)$$

where

$$\Delta T_k(s) = T_k(s) - T_c(s) \quad \text{for } k = 1, 2, 3. \quad (15)$$
<table><caption>TABLE I: Parameter values for Havriliak-Negami models $M_{12}^{HN}(s)$ and $M_{23}^{HN}(s)$.</caption><thead><tr><th></th><th></th><th>ν<sub>1</sub></th><th>ν<sub>2</sub></th><th>λ</th><th>ν<sub>0</sub></th><th>K</th><th>J<sub>ident</sub></th><th>J<sub>val</sub></th></tr></thead><tbody><tr><td rowspan="2">M<sub>12</sub><sup>HN</sup>(s)</td><td>pre ECC</td><td>1.50</td><td>-0.32</td><td>3.24</td><td>0.24</td><td>0.87</td><td>7.6 10<sup>-4</sup></td><td>7.3 10<sup>-4</sup></td></tr><tr><td>ECC</td><td>1.65</td><td>1.20</td><td>3.40</td><td>-10<sup>-3</sup></td><td>0.67</td><td>2.3 10<sup>-3</sup></td><td>3.5 10<sup>-3</sup></td></tr><tr><td rowspan="2">M<sub>23</sub><sup>HN</sup>(s)</td><td>pre ECC</td><td>1.48</td><td>1.31</td><td>8.79</td><td>0.11</td><td>2.27</td><td>9.5 10<sup>-4</sup></td><td>8.8 10<sup>-4</sup></td></tr><tr><td>ECC</td><td>1.92</td><td>-0.06</td><td>1.78</td><td>0.09</td><td>1.25</td><td>7.2 10<sup>-3</sup></td><td>12 10<sup>-3</sup></td></tr></tbody></table>

<table><caption>TABLE II: Parameter values for rational models $M_{12}^{rat}(s)$ and $M_{23}^{rat}(s)$.</caption><thead><tr><th></th><th></th><th>a<sub>0</sub></th><th>a<sub>1</sub></th><th>b<sub>1</sub></th><th>b<sub>2</sub></th><th>b<sub>3</sub></th><th>J<sub>ident</sub></th><th>J<sub>val</sub></th></tr></thead><tbody><tr><td rowspan="2">M<sub>12</sub><sup>rat</sup>(s)</td><td>pre ECC</td><td>0.03</td><td>1.62</td><td>1.50</td><td>0.25</td><td>0.03</td><td>6.0 10<sup>-2</sup></td><td>6.0 10<sup>-2</sup></td></tr><tr><td>ECC</td><td>-49.86</td><td>2923</td><td>9995</td><td>2.2</td><td>1017</td><td>5.1 10<sup>-2</sup></td><td>7.6 10<sup>-2</sup></td></tr><tr><td rowspan="2">M<sub>23</sub><sup>rat</sup>(s)</td><td>pre ECC</td><td>0.14</td><td>8.00</td><td>3.57</td><td>0.02</td><td>0.02</td><td>5.5 10<sup>-2</sup></td><td>4.9 10<sup>-2</sup></td></tr><tr><td>ECC</td><td>-276.5</td><td>4353</td><td>1000</td><td>3.718</td><td>650.7</td><td>1.3 10<sup>-1</sup></td><td>1.9 10<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 6: Validation data and modeling error of $H_{12}(s)$ before ECC</figcaption></figure>

A. *Comparison between Havriliak-Negami models and rational models*

For system identification of the transfers (13) and (14), models are identified using a criterion of output error type by using the nonlinear simplex algorithm [23, 24]. Two different models are used: one using the Havriliak-Negami function

$$M_j^{HN}(s) = K s^{\nu_0} (s^{\nu_1} + \lambda)^{\nu_2} \quad (16)$$

where the parameter vector $[\nu_1 \nu_2 \lambda \nu_0 K]$ is optimized and one using a rational transfer function

$$M_j^{rat}(s) = \frac{A(s)}{B(s)} = \frac{a_0 + a_1 s}{1 + b_1 s + b_2 s^2 + b_3 s^3} \quad (17)$$

where the parameter vector $[a_0 a_1 b_1 b_2 b_3]$ is optimized. In models (16) and (17), $j$ stands for 12 or 23; therefore, $M_{12}^{HN}(s)$ and $M_{12}^{rat}(s)$ model the dynamics of $H_{12}(s)$, while $M_{23}^{HN}(s)$ and $M_{23}^{rat}(s)$ model the dynamics of $H_{23}(s)$. Measures are collected at sampling frequency $F_s = 50$ Hz, therefore the Havriliak-Negami function is then approximated on the frequency band $[\omega_b, \omega_h] = [1.57, 157]$ (rad/s).

The estimated parameters are presented on Table I for Havriliak-Negami transfer functions and on Table II for rational transfer functions. Note that the collected data are divided into two parts, the first part being used for the model parameter estimation (see the column of the criterion evaluation $J_{ident}$) and the other part being used for the validation purpose (see the
<table><caption>TABLE III: Parameter values for fractional models $M_{12}^{frac}(s)$ and $M_{23}^{frac}(s)$.</caption><thead><tr><th></th><th></th><th>a<sub>0</sub></th><th>a<sub>1</sub></th><th>b<sub>1</sub></th><th>b<sub>2</sub></th><th>&nu;</th><th>J<sub>ident</sub></th><th>J<sub>val</sub></th></tr></thead><tbody><tr><td rowspan="2">M<sub>12</sub><sup>frac</sup>(s)</td><td>pre ECC</td><td>0.46</td><td>0.36</td><td>-0.23</td><td>0.19</td><td>0.34</td><td>8.3 10<sup>-4</sup></td><td>8.1 10<sup>-4</sup></td></tr><tr><td>ECC</td><td>0.51</td><td>0.05</td><td>-0.14</td><td>0.017</td><td>0.65</td><td>8.0 10<sup>-3</sup></td><td>16.0 10<sup>-3</sup></td></tr><tr><td rowspan="2">M<sub>23</sub><sup>frac</sup>(s)</td><td>pre ECC</td><td>1.17</td><td>0.75</td><td>-0.24</td><td>0.15</td><td>0.07</td><td>9.5 10<sup>-4</sup></td><td>8.9 10<sup>-4</sup></td></tr><tr><td>ECC</td><td>0.97</td><td>0.25</td><td>-0.23</td><td>0.20</td><td>0.52</td><td>83.0 10<sup>-3</sup></td><td>14.0 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 7: Validation data and modeling error of $H_{12}(s)$ during ECC</figcaption></figure>

criterion evaluation $J_{val}$). Note that the precision of the rational and fractional models are compared by using the optimization criterion $J$ ($J_{ident}$ for identification data and $J_{val}$ for validation data) defined as:

$$J_{ident/val} = \frac{1}{N} \sum_{k=1}^{N} (Y(k) - \hat{Y}(k))^2 \quad (18)$$

where $N$ denotes the number of samples of the real output response $Y$, $\hat{Y}$ being the estimated response. Note that the criterion values are always minimal in the case of fractional models (100 times lower). With identical number of parameters, the precision of fractional models is therefore better than the one of rational models. It has been shown in [20] that to get the same level of precision, a rational model needs to have at least two times more parameters than a fractional model.

Figures 6, 7, 8 and 9 show the results of the estimated models on validation data. The blue curves represents the exciting input signal ($\Delta T_k$), the black ones the system response ($\Delta T_{k+1}$), the red curves the model response estimated with Havriliak-Negami function, and the green ones those of the model estimated with a rational transfer function. These figures confirm the results presented on Tables I and II when comparing the criterion values ($J_{val}$) and also show a better precision of Havriliak-Negami models versus rational models.

## B. Comparison between Havriliak-Negami models and fractional order ones

In order to highlight the interest of Havriliak-Negami function for system identification of fractional models, transfers $H_{12}(s)$ and $H_{23}(s)$ are now estimated by using more classic fractional transfer functions of explicit form such as:

$$M_j^{frac}(s) = \frac{a_0 + a_1 s^{\nu}}{1 + b_1 s^{\nu} + b_2 s^{2\nu}} \quad (19)$$
<figure><img src="image_10.png" /><figcaption>Fig. 8: Validation data and modeling error of $H_{23}(s)$ before ECC</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9: Validation data and modeling error of $H_{23}(s)$ during ECC</figcaption></figure>

with $j = 12$ or $23$ and $\nu$ being the commensurate order<fnref n="2" />. $M_{12}^{frac}(s)$ models the dynamics of $H_{12}(s)$, while $M_{23}^{frac}(s)$ models the dynamics of $H_{23}(s)$.

<footnote n="2">²When all differentiation orders are integer multiple of a single order $\nu$, called commensurate, the searching domain of the commensurate order is held for $\nu \in (0, 2)$, corresponding to the stability domain (see e.g. Matignon's theorem [25])</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 10: Validation data and modeling error of $H_{12}(s)$ before ECC</figcaption></figure>

The parameter vector $[a_0 a_1 b_1 b_2 \nu]$ is estimated using the output error function $oe$ provided by the CRONE toolbox developed in Matlab® (see [26] for the latest developments on system identification of fractional models using the CRONE toolbox). The $oe$ function is implemented using the gradient based optimization algorithm of Levenberg-Marquardt. The estimated parameters are presented on Table III.

As the models are identified as black-box models, it is difficult to give a physical or physiological interpretation of these parameters. Nevertheless, it is interesting to note that models $M_j^{frac}(s)$ have an estimated commensurate order around 0.5 for measures carried out during the ECC whatever the transfer $H_{12}(s)$ or $H_{23}(s)$.

Figures 10, 11, 12 and 13 show the results of the estimated models on validation data. The blue curves represents the exciting input signal ($\Delta T_k$), the black ones the system response ($\Delta T_{k+1}$), the red curves the model response estimated with Havriliak-Negami function, and the green ones those of the model estimated with a rational transfer function. These figures confirm the results presented on Tables I and III when comparing the criterion values ($J_{val}$). Havriliak-Negami models have better precision compared to classic fractional models with the same number of parameters. Releasing the commensurate order constraint would help the precision however the number of parameters will double. These figures and the criterion value table show that the obtained results with Havriliak-Negami models $M_j^{HN}(s)$ are better than those obtained with fractional models $M_j^{frac}$ for measures realized before ECC, that is to say, in natural operating condition of the lungs. As for the measures achieved during the ECC, the results in terms of precision are equivalent.

### C. Comparison of the models in frequency domain

The results obtained by comparing the time performances of the proposed models are to be moderated with the frequency study of the input signals. Indeed, the power spectral density of the input signal $\Delta T_1$ and $\Delta T_2$ show that most of its energy is confined around the frequency $0.2 \, \text{Hz}$ or $1.89 \, \text{rad/s}$, the operating conditions being natural or during an ECC. This operating frequency does explain the quasi-periodic form of the input signals (see Figure 14), a form that is due to the change of air ventilation cycles (filling/inflation) inside the lungs. The identified models are therefore valid only around this frequency.

Furthermore, in order to complete the frequency analysis on the results, the Bode diagrams of the models have been studied and plot on Figures 15, 16, 17 and 18.

System identification works well when the input signals are persistently exciting. Also, the power spectral density of the input signals reveal the exciting frequency of the input signals are around $1.9 \, \text{rad/s}$. Consequently, the system identification results are bounded to the exciting input. Enriching the input signal to be applied to the lungs to other frequencies would be better, but might damage the lungs. Nevertheless, it could be studied in future works. The Bode diagrams do confirm that the different models are indeed close around the exciting frequency of $1.9 \, \text{rad/s}$ (or $0.2 \, \text{Hz}$). Moreover, some of the models
<figure><img src="image_13.png" /><figcaption>Fig. 11: Validation data and modeling error of $H_{12}(s)$ during ECC</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12: Validation data and modeling error of $H_{23}(s)$ before ECC</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 13: Validation data and modeling error of $H_{23}(s)$ during ECC</figcaption></figure>

appear to have differentiator behavior in the frequency band of interest. This observation calls out the input/output causality principle of the studied system and of transfers $H_{12}(s)$ and $H_{23}(s)$ as defined at the beginning of section IV.

The system causality is a perspective to study and is not completely known. Actually, during this experience, first the lung is insufflated with fresh air, then this same air is sucked after being reheated with the contact of warmer tissues of the lung. The way of the air in movement thus varies during the experience in function of the insufflation/suction cycles. The air moves inside the lungs during the insufflation and then moves outside while the suction. A conceivable hypothesis to explain the system causality could be to consider that it varies according to the direction of the air in the lung.

V. CONCLUSION

Eight models of bronchus have been identified by considering two different states of the system with or without extracorporeal circulation. The estimated fractional models, including the Havriliak-Negami models as they also are fractional models of implicit form, enable rendering in an accurate manner the temperature variations inside the lungs.

The comparison between the models using the Havriliak-Negami function and the rational ones has shown that with a same number of parameters, the fractional models are more accurate to model the studied system. This observation confirms the hypothesis that the system dynamics are governed by fractional order differentiation equations.

The comparison of fractional order model performance has shown that, with identical number of parameters, Havriliak-Negami models are more accurate in natural operating conditions, and identically accurate during ECC. These observations do reinforce the interest of using Havriliak-Negami functions for black-box system identification.

Moreover, a significant variation of model parameters has been observed between both system operating conditions (before and during ECC). For the moment, these parameter variations have not been interpreted physically nor physiologically. Nevertheless, the parameter variations are important in the case of classic fractional transfer functions and the commensurate order tends towards 0.5 for measures carried out during ECC which also confirms fractional differential equations as the solving of the heat equation does reveal a fractional order of 0.5.

The perspectives of this work are on the investigation of physical or physiological sense of the parameter variations. The possibility to use an analytical model, such as an arborescent arrangement or a cascade arrangement as described in section II, should be studied in order to give an easier physical interpretation. Finally, a study on the system causality should also be considered in order to better understand the system dynamics and to better identify it.

REFERENCES

[1] A. Oustaloup, *La dérivation non-entière : théorie, synthèse et applications.* Paris: Hermès, 1995.
<figure><img src="image_16.png" /><figcaption>Fig. 14: Power spectral density of a) input $\Delta T_1$ (pre ECC), b) input $\Delta T_1$ (during ECC), c) input $\Delta T_2$ (pre ECC), and d) input $\Delta T_2$ (during ECC)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 15: Bode diagram of the estimated models to model *H*₁₂, before ECC</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16: Bode diagram of the estimated models to model *H*₁₂, during ECC</figcaption></figure>

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

<footnote>[3] A. Maachou, R. Malti, P. Melchior, J. Battaglia, A. Oustaloup, and B. Hay, "Nonlinear thermal system identification using fractional Volterra series," *Control Engineering Practice*, vol. 29, pp. 50-60, 2014.</footnote>
<figure><img src="image_19.png" /><figcaption>Fig. 17: Bode diagram of the estimated models to model *H*₂₃, before ECC</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18: Bode diagram of the estimated models to model *H*₂₃, during ECC</figcaption></figure>

<footnote>[4] S. Rodrigues, N. Munichandraiah, and A.-K. Shukla, “A review of state of charge indication of batteries by means of A.C. impedance measurements,” *Journal of Power Sources*, vol. 87, no. 1-2, pp. 12–20, 2000.</footnote>

<footnote>[5] N. Heymans and J. Bauwens, “Fractal rheological models and fractional differential equations for viscoelastic behavior,” *Rheologica Acta*, vol. 33, p. 219, 1994.</footnote>
[6] A. Benchellal, S. Bachir, T. Poinot, and J.-C. Trigeassou, "Identification of a non-integer model of induction machines," in *1st IFAC Workshop on Fractional Differentiation and its Applications (FDA)*, Bordeaux, France, 2004, pp. 400-407.

[7] B. Mandelbrot, *Les objets fractals, forme, hasard et dimension*. Flammarion, 1975.

[8] C. M. Ionescu and R. D. Keyser, "Relations between fractional-order model parameters and lung pathology in chronic obstructive pulmonary disease," *IEEE Transactions on Biomedical Engineering*, vol. 56, no. 4, pp. 978-987, 2009.

[9] C. Ionescu, J. Machado, and R. De Keyser, "Modeling of the lung impedance using a fractional-order ladder network with constant phase elements," *IEEE Transactions on Biomedical Circuits and Systems*, vol. 5, no. 1, pp. 83-89, 2011.

[10] R. L. Magin, *Fractional Calculus in Bioengineering*. Begell House Publisher, Inc. Connecticut, 2006.

[11] R. L. Magin and M. Ovadia, "Modeling the cardiac tissue electrode interface using fractional calculus," *Journal of Vibration and Control*, vol. 19, pp. 1431-1442, 2009.

[12] A. Oustaloup, *Diversity and Non-integer Differentiation for System Dynamics*. Wiley-ISTE, 2014.

[13] P. Melchior, M. Pellet, Y. Abdelmoumen, and A. Oustaloup, "System identification of thermal transfers inside the lungs using fractional models," in *7th Vienna International Conference on Mathematical Modelling*, vol. 45, no. 2, 2012, pp. 571-576, iFAC Proceedings Volumes.

[14] Z. Dai, Y. Peng, H. A. Mansy, R. H. Sandler, and T. J. Royston, "A model of lung parenchyma stress relaxation using fractional viscoelasticity," *Medical Engineering and Physics*, vol. 37, no. 8, pp. 752-758, 2015.

[15] C. Ionescu, R. Hodrea, and R. De Keyser, "Variable time-delay estimation for anesthesia control during intensive care," *IEEE Transactions on Biomedical Engineering*, vol. 58, no. 2, pp. 363-369, 2011.

[16] L. Sommacal, "Synthèse de la fonction d'Havriliak-Negami pour l'identification temporelle par modèle non entier et modélisation du système musculaire," Ph.D. dissertation, Université Bordeaux 1, Talence, 2007.

[17] M. Pellet, P. Melchior, Y. Abdelmoumen, and A. Oustaloup, "Fractional thermal model of the lungs using havriliak-negami function," in *International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, 2011 ASME/IEEE International Conference on Mechatronic and Embedded Systems and Applications, Parts A and B*, vol. 3, no. 54808, 2011, pp. 805-810.

[18] M. Pellet, "Caractérisation non entière de systèmes biologiques : application au muscle squelettique et au poumon," Ph.D. dissertation, Université Bordeaux 1, Talence, France, July 2013.

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

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

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

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

[23] M. Subrahmanyam, "An extension to the simplex method to constrained nonlinear optimization," *Journal of Optimization Theory and Applications*, vol. 62, no. 2, pp. 311-319, Aug 1989.

[24] D. Woods, "An interactive approach for solving multi-objective optimization problems (interactive computer, nelder-mead simplex algorithm, graphics)," Ph.D. dissertation, Rice University, 1985.

[25] D. Matignon and B. d'Andréa-Novel, "Some results on controllability and observability of finite-dimensional fractional differential systems," in *IMACS*, vol. 2. Lille, France: IEEE-SMC, July 1996, pp. 952-956.

[26] R. Malti and S. Victor, "Crone toolbox for system identification using fractional order models," in *17th IFAC Symposium on System Identification (SYSID 2015)*, vol. 48, no. 28, Beijing, China, October 2015, pp. 769-774.

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

**Stéphane Victor** was born in Germany in 1983. He graduated from the Bordeaux INP/ENSEIRB-MATMECA and Ecole Polytechnique de Montrã©al engineering schools in 2006, and obtained his M.Sc. and Ph.D. degrees in Automatic control from Universitã© de Bordeaux in 2006 and 2010. 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 the area of fractional differentiation and its applications in automatic control, system identification, thermal systems, autonomous vehicles and motion planning with flatness.
<figure><img src="image_22.png" /></figure>

**Pierre Melchior** was born in France on 1961. He graduated from the University of Bordeaux, Talence, France, and received the Ph.D. degree in control engineering from Paul Sabatier University, Toulouse, France, in 1989 and the Accreditation to supervise research (HDR) from the University of Bordeaux, in 2014. He is currently Professor of control theory and a Researcher with Laboratory IMS, University of Bordeaux, Talence. He is head of the CRONE team of IMS Automatic Control Group. His research interests include systems, control theory, robotics, and biomechatronics (e.g., fractional systems, fractional path planning, fractional motion planning, and muscle modeling).

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

**Mathieu Pellet** was born in France in 1985. He graduated from université Bordeaux I in 2008, and received his Ph.D. degree in 2013 from université Bordeaux I. His research interests are in the area of fractional differentiation and its applications in automatic control, biological systems, system identification, thermal systems.

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

**Alain Oustaloup** was born in France in 1950. He received his Engineering diploma from ENSEIRB in 1973, and his Doctor of Engineering and Doctor of Science in 1975 and 1981 from the University of Bordeaux. He is currently Professor at Bordeaux INP/ENSEIRB-MATMECA. He founded the IMS-CRONE team at the University of Bordeaux. He is currently a Research Evaluator at the French Ministry of Education and Research. His research interests include fractional differentiation, its synthesis, and its applications in Engineering sciences, particularly in automatic control. He conceived CRONE control (French acronym for "Commande Robuste d'Ordre Non Entier"). He is author and co-author of five books including the reference Oustaloup (1995). He was awarded the Afcet Trophy in 1995, the CNRS Silver Medal in 1997, and the Grand Prix Lazare Carnot in 2011 from the french Académie des Sciences.