# Time-Domain Fractional Behaviour Modelling with Rational Non-Singular Kernels

Jocelyn Sabatier, Christophe Farges

► To cite this version:

Jocelyn Sabatier, Christophe Farges. Time-Domain Fractional Behaviour Modelling with Rational Non-Singular Kernels. Axioms, 2024, 13 (2), pp.99. 10.3390/axioms13020099. hal-05455108

HAL Id: hal-05455108

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

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
# Time-Domain Fractional Behaviour Modelling with Rational Non-Singular Kernels

Jocelyn Sabatier * and Christophe Farges

IMS Laboratory, Bordeaux University, UMR 5218 CNRS, 351 Cours de la Libération, 33405 Talence, France; christophe.farges@u-bordeaux.fr

* Correspondence: jocelyn.sabatier@u-bordeaux.fr

**Abstract:** This paper proposes a solution to model fractional behaviours with a convolution model involving non-singular kernels and without using fractional calculus. The non-singular kernels considered are rational functions of time. The interest of this class of kernel is demonstrated with a pure power law function that can be approximated in the time domain by a rational function whose pole and zeros are interlaced and linked by geometric laws. The Laplace transform and frequency response of this class of kernel is given and compared with an approximation found in the literature. The comparison reveals less phase oscillation with the solution proposed by the authors. A parameter estimation method is finally proposed to obtain the rational kernel model for general fractional behaviour. An application performed with this estimation method demonstrates the interest in non-singular rational kernels to model fractional behaviours. Another interest is the physical interpretation fractional behaviours that can be implemented with delay distributions.

**Keywords:** fractional behaviours; convolution model; impulse response; non-singular rational kernels

**MSC:** 45E10; 26A33

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

check for
updates

**Citation:** Sabatier, J.; Farges, C. Time-Domain Fractional Behaviour Modelling with Rational Non-Singular Kernels. *Axioms* 2024, 13, 99. https://doi.org/10.3390/axioms13020099

Academic Editor: Trushit Patel

Received: 15 December 2023

Revised: 23 January 2024

Accepted: 26 January 2024

Published: 31 January 2024

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

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Integral operators with power-law-based kernels are usually used to model fractional behaviour of a physical origin [1–3]. This is the case with the fractional order integration operator defined with the Riemann–Liouville definition [4]

$$ 
\begin{aligned}
I^{\gamma}[u(t)] &= \frac{1}{\Gamma(\gamma)} \int_{0}^{t} \frac{1}{(t-\tau)^{1-\gamma}} u(\tau) d\tau = \int_{0}^{t} \eta_{PL}(t-\tau) u(\tau) d\tau \\
&= \int_{0}^{t} \eta_{PL}(t) u(\tau-\tau) d\tau = \eta_{PL}(t) * u(t)
\end{aligned} 
\quad (1) $$

which involves a power law kernel $\eta_{PL}(t) = t^{\gamma-1}/\Gamma(\gamma)$. This power kernel can also be viewed as the impulse response or the Green function of the fractional integration operator. This operator is then used to define more complex models such as pseudo-state space descriptions or fractional differential equations or often used to capture fractional behaviours [5–7]. In a general way, we will say that a system exhibits a fractional behaviour if its time or frequency response can be fitted by a power law behaviour in a given time or frequency range. More precisely, as shown in [3], this power law behaviour can be highlighted using spectral density and autocorrelation functions of the system input-output signals.

However, the above cited kernel is singular in the sense that it tends towards infinity as the time tends towards 0. As explained in [8], the definition of fractional operators with singular kernels leads to cumbersome mathematical expressions and the consequent complications in the solutions of the fractional order differential equations [8]. This class of kernel is also the origin of several limitations associated with fractional-integration-based models that have been described in the literature (initialisation issues, infinitely
small and high time constants, infinite memory length, etc.) [9–13]. It is for these reasons that some authors have suggested definitions of fractional operators using non-singular kernels [14–17], as this class of kernel permits most of the above-mentioned limitations to be solved.

But when it is realized that fractional behaviours (those produced by physical phenomena) and fractional models are two different concepts, since the former designates a property of a physical system and the latter designates a model class, among a set of model classes that capture fractional behaviours, then a whole new field opens up, covering (among other things) the search for new models intended for fractional behaviours, the study of these new models’ properties, and the search for new identification methods. Some studies conducted from this perspective have already been reported in the literature [18,19]. The present paper pursues this line of inquiry.

In this paper, we propose to model fractional behaviour without using fractional calculus and without resorting to a frequency representation of the models and of their approximation. Classically and in practice, modelling fractional behaviour using a fractional model identified from time data requires the following steps:

- Select a fractional model represented by a transfer function or a pseudo state space description with a given complexity (and thus a set of associated parameters);

- Approximate the fractional operators in this model to obtain an equivalent integer system or a truncated series that permits time simulations, which often leads to a very high-dimensional integer model;

- Implement a minimization method to find the model parameters that reduce the error between the response model response and the data.

To use the resulting optimal model (for simulation, control, etc.), the high dimensional approximation previously mentioned is again used.

In this paper, we model fractional behaviours using integral operators with rational kernels, which will have the following two advantages:

- In a modelling approach, the approximation step previously mentioned is avoided and replaced by the discretization of Equation (1), which is simplified as the proposed rational kernel is not singular, contrary to $\eta_{PL}(t)$ in Equation (1);

- In comparison to degree $n$ polynomials, the approximation by a degree $n$ rational function produces a very small approximation error.

As an example of the latter claim, one can mention Newman’s theorem [20], which shows that the function $|t|$ can be approximated on $[-1, 1]$ by rational functions of degree $n$ with an error $O(e^{-K\sqrt{n}})$, whereas it is known that an error larger than $K/n$ takes place with degree $n$ polynomials.

Moreover, for the power law function $t^{-\nu}$, it is demonstrated in the paper that a particular distribution of zeros and poles in the rational kernel allows an accurate fitting on a time domain range, with a very low number of parameters and fewer frequency oscillations in comparison to a distribution of exponential functions as conventionally used in the literature [21–27].

The paper is thus organized as follows. Section 2 deals with the approximation of the power law kernel $Kt^{-\nu}$ by a rational kernel with a particular distribution of poles and zeros. The advantages of such an approximation are highlighted, which justifies their use to model fractional behaviours. Section 3 gives the Laplace transform and frequency response of these kernels. Some improvements in the kernel definitions are also proposed. In Section 4, the proposed kernel class is used to model a fractional behaviour. The goal of this section is not to propose a robust identification methodology, but to show the efficiency of the proposed class of kernel using only a nonlinear optimization algorithm. In Section 5, we show that the rational kernels considered can also be represented by a distribution of delays, which is close to physical phenomena such as diffusion, adsorption or aggregation, which are known to produce power law (or fractional) behaviours. For these phenomena,
the delays in question are the times required by the reacting agents (atoms, molecules, people, etc.) to diffuse up to a given abscissa, to absorb or even to aggregate.

## 2. Power Law Kernel Approximation with a Rational Kernel

In this section, we consider the approximation of a pure power law function of the form

$$ \eta_v(t) = Kt^{-\nu} \quad (2) $$

Applying a decimal logarithm to the function $\eta_v(t)$ leads to:

$$ \log(\eta_v(t)) = \log(K) - \nu \log(t). \quad (3) $$

We propose to approximate $\eta_v(t)$ by the rational kernel

$$ \eta(t) = C_0 \frac{\prod_{j=1}^{N} \left( \frac{t}{t'_j} + 1 \right)}{\prod_{j=1}^{N} \left( \frac{t}{t'_j} + 1 \right)} \quad t'_j > 0, \quad t_j > 0. \quad (4) $$

The decimal logarithm of the function $\eta(t)$ is

$$ \log(\eta(t)) = \log(C_0) + \sum_{j=1}^{N} \log\left(\frac{t}{t'_j} + 1\right) - \sum_{j=1}^{N} \log\left(\frac{t}{t'_j} + 1\right) \quad (5) $$

At $t = 1$, according to (2),

$$ \eta_v(1) = K \quad (6) $$

and, according to (4),

$$ \eta(1) = C_0 \frac{\prod_{j=1}^{N} \left( \frac{1}{t'_j} + 1 \right)}{\prod_{j=1}^{N} \left( \frac{1}{t'_j} + 1 \right)}. \quad (7) $$

Consequently, imposing $\eta(1) = \eta_v(1)$ leads to

$$ C_0 = K \frac{\prod_{j=1}^{N} \left( \frac{1}{t'_j} + 1 \right)}{\prod_{j=1}^{N} \left( \frac{1}{t'_j} + 1 \right)}. \quad (8) $$

Figure 1 represents the terms $\eta_z(t) = \log(\frac{t}{t'_j} + 1)$ and $\eta_p(t) = -\log(\frac{t}{t'_j} + 1)$ from Equation (5) for $t'_j = t_j = 1$ as a function of $\log(t)$. This figure shows that for $t > (t'_j = 1)$ and $t > (t_j = 1)$, the functions $\eta_z(t)$ and $\eta_p(t)$, respectively, tend towards lines of slope of 1 and -1. Based on the work of Manabe [21], Oustaloup [24] and Charef [26], the function defined by Equation (3), whose slope in logarithm is $-\nu$, can thus be approximated by an alternation of functions $\eta_p(t)$ and $\eta_z(t)$ and by imposing the appropriate values of $t'_j$ and $t_j$.

Figure 2 shows how the functions $\eta_z(t)$ and $\eta_p(t)$ are interlaced to obtain the desired slope, assuming that it is required to approximate Equation (2) on the time interval $[t_l, t_h]$. It is moreover assumed that the times $t'_j$ and $t_j$, which enter into the definitions of $\eta_z(t)$ and $\eta_p(t)$, meet the following equalities:

$$ \log(\alpha) = \log(t'_j) - \log(t_j) \quad \text{and thus} \quad \alpha = \frac{t'_j}{t_j} \quad (9) $$
$$ \log(\eta) = \log(t_{j+1}) - \log(t'_j) \quad \text{and thus} \quad \beta = \frac{t_{j+1}}{t'_j} \qquad (10) $$

$$ \log(\alpha) + \log(\beta) = \log(t_{j+1}) - \log(t_j) = \log(t'_{j+1}) - \log(t'_j) \quad \text{and thus} \quad r = \alpha\beta = \frac{t_{j+1}}{t_j} = \frac{t'_{j+1}}{t'_j} \qquad (11) $$

<figure><img src="image_4.png" /><figcaption>Figure 1. Functions $\eta_z(t) = \log(t/t'_j + 1)$ and $\eta_p(t) = -\log(t/t'_j + 1)$ for $t'_j = t_j = 1$, as functions of $\log(t)$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2. Interlacing of function $\eta_p(t)$ and $\eta_z(t)$ (blue line) used in the algorithm for the approximation of an affine function of slope $v$ (red line).</figcaption></figure>
If an amount of N functions $\eta_z(t)$ and $\eta_p(t)$ are used in the approximation, then using Figure 2, the following equality holds:

$$N(\log(\alpha) + \log(\beta)) = \log(t_h) - \log(t_l) \text{ and thus } (\alpha\beta)^N = r^N = \frac{t_h}{t_l} \quad (12)$$

So that function (2) and its approximation have the same slope, the values of times $t_j$ and $t'_j$ must be meet the following equalities:

$$\frac{A}{\log(\alpha) + \log(\beta)} = -\nu \text{ and } \frac{A}{\log(\alpha)} = -1 \quad (13)$$

that can be combined into:

$$\frac{\log(\alpha)}{\log(\alpha) + \log(\beta)} = \nu \text{ thus } \alpha = (\alpha\beta)^{\nu}. \quad (14)$$

The efficiency of the above approximation method is illustrated by Figures 3 and 4 for the power law function (2) with $\nu = 0.5$ and $K = 10$. The following parameters were also used in this application, $N = 7$, $t_l = 10^{-3}$, $t_h = 10^4$, and the following steps are followed:

1. Compute $r = \sqrt[N]{\frac{t_h}{t_l}} = 10$;

2. Compute $\alpha = r^{\nu} = 3.162$ and $\beta = \frac{r}{\alpha} = 3.162$;

3. Compute $t'_1 = t_l \times \sqrt{\beta} = 5.6234 \cdot 10^{-3}$ and the other $t'_j$ and $t_j$ using Equations (9) and (10).

<figure><img src="image_6.png" /><figcaption>Figure 3. Comparison of $\log(\eta(t))$ and $\log(\eta_\nu(t))$ with a logarithmic time scale.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4. Comparison of $\eta_v(t)$ and $\eta(t)$.</figcaption></figure>

As in the works of Manabe [21], Oustaloup [24] and Charef [26], the kernel and the recursive relation we propose thus permit a power law behaviour to be fitted in a given time range, where [21,24,26] fit it on a frequency range. The interest in working directly in the time domain will appear in the modelling application of Section 4, as usually data produced by a real system to model are time data.

The analysis carried out in this section permits us to conclude that the new kernel (4) we propose can be viewed as a serious candidate for extensions of fractional operator definition as those in [28–31]. Perhaps the operators that will be formed with this class of kernels will not be able to be qualified as fractional derivation or integration operators in the strict sense. But, as stated in the introduction, the primary objective is to have operators capable of realistically taking into account the memory of real fractional behaviours without falling under the cost of the limitations mentioned in the introduction.

### 3. Laplace Transform and Frequency Response of the Obtained Rational Kernel

A frequency representation of kernel (4) can be obtained after its partial fraction decomposition, which can be written as follows:

$$ \eta(t) = C'_{0} + \sum_{j=1}^{N} \frac{a_{j}}{\left(\frac{t}{t_{j}} + 1\right)} = C'_{0} + \sum_{j=1}^{N} \frac{a_{j}t_{j}}{(t + t_{j})}. \quad (15) $$
Let $F(s)$ be the Laplace transform of the time function $f(t)$, namely $F(s) = \int_{0}^{\infty} f(\tau)e^{-s\tau} d\tau = \mathcal{L}\{f(t)\}$, with $s$ being the Laplace variable. As [32,33]:

$$ \mathcal{L}\left\{\frac{1}{t+b}\right\} = e^{sb}Ei_1.sb \quad \text{with} \quad Ei_1(x) = \int_x^\infty \frac{e^{-\zeta}}{\zeta}d\zeta \qquad (16) $$

then, by analogy,

$$ \mathcal{L}\{\eta(t)\} = \eta(s) = \frac{C_0'}{s} + \sum_{j=1}^{N} a_j t_j e^{st_j} Ei_1(st_j). \qquad (17) $$

Note that for $\mathcal{Re}(x) > 0$ [33], so we also have:

$$ Ei_1(x) = \int_1^\infty \frac{e^{-\zeta x}}{\zeta} d\zeta = \int_0^1 \frac{e^{-\zeta/x}}{\zeta} d\zeta. \qquad (18) $$

See Appendix A for the evaluation of $Ei_1(s)$, if $s = j\omega$ (frequency response). Figure 5 shows the Bode diagram of $\eta(j\omega)$ for the same example as in Section 2.

<figure><img src="image_8.png" /><figcaption>Figure 5. Bode diagram of Laplace transform of function $\eta(t)$.</figcaption></figure>

A comparison is now proposed between the approximation $\eta(t)$ given by Equation (4) and the well-known [21–27] approximation based on a recursive distribution of the exponential that is defined in the time domain by:

$$ \chi(t) = K_0 + \sum_{j=1}^{N} a_j e^{-\frac{t}{r_j \omega_1}}. \qquad (19) $$

The inverse Laplace transform of the relation of Equation (19) is given by:

$$ \chi(s) = \frac{K_0 \prod_{j=1}^{N} \left(1 + \frac{s}{\omega_j'}\right)}{s \prod_{j=1}^{N} \left(1 + \frac{s}{\omega_j}\right)} \qquad (20) $$

with

$$ \omega'_{j+1} = r\omega'_j \quad \omega_j = r^\nu \omega'_j. \qquad (21) $$

This comparison is shown in Figure 6, in which parameters $\nu$ and $N$ are, respectively, chosen equal to 0.3 and 7 for Equations (17) and (20). For a fair comparison, the two
approximations must cover the same frequency range. The parameters $K_0$, $r$ and $\omega'_1$ in Equation (20) and $C_0$, $r$ and $t_1$ in Equation (17) are thus defined by:

$$C'_0 = 0.773, t_1 = 2.2387 * 10^{-3} \text{ s, } r = 10 \text{ for Equation (17)},$$

$$K_0 = 0.522, \omega_1 = 1.1115 * 10^{-4} \text{ rd/s, } r = 9.8 \text{ for Equation (20).}$$

A comparison with the fractional order transfer function

$$H_f(s) = \frac{H_0}{s} \left(1 + \frac{s}{\omega_l}\right)^{\nu} \left(1 + \frac{s}{\omega_h}\right)^{\nu}. \qquad (22)$$

with $H_0 = 0.5$, $\omega_l = 45 * 10^{-6} \text{ rd/s}$, $\omega_h = 550 \text{ rd/s}$ is also shown in Figure 6.

<figure><img src="image_9.png" /><figcaption>Figure 6. Comparison of the rational approximation (17), the exponential approximation (20) and the transfer function (22).</figcaption></figure>

This figure shows that:

- The gain diagrams are similar for the three functions;

- The rational function leads to slower phase changes, mainly in high frequencies (but the impact on the time response is low as the gain is very small in this frequency range);

- The phase oscillations are much less pronounced with the rational function than with the exponential sum.

Other fractional behaviours can be obtained by adding an extra pole $t_{N+1} = rt_N$ in Equation (4). The corresponding Bode diagram is given in Figure 7 and there exhibits a constant behaviour in low frequency.

The behaviour at the low corner frequencies can be made more abrupt by adding an exponential function in Equation (4), which becomes:

$$\eta(t) = C_0 \frac{\prod_{j=1}^{N} \left( \frac{t}{t'_j} + 1 \right)}{\prod_{j=1}^{N} \left( \frac{t}{t'_j} + 1 \right)} e^{-at}. \qquad (23)$$

The decimal logarithm of function $\eta(t)$ is defined by

$$\log(\eta(t)) = \log(C_0) - \frac{at}{\ln(10)} + \sum_{j=1}^{N} \log\left(\frac{t}{t'_j} + 1\right) - \sum_{j=1}^{N} \log\left(\frac{t}{t'_j} + 1\right). \qquad (24)$$
Partial fraction decomposition of Equation (10) leads to:

$$ \eta(t) = \left( C'_{0} + \sum_{j=1}^{N} \frac{a_j}{\left( \frac{t}{t_j} + 1 \right)} \right) e^{-at}. \qquad (25) $$

As

$$ \mathcal{L} \left\{ \frac{e^{-at}}{\left( \frac{t}{t_j} + 1 \right)} \right\} = t_j e^{t_j(a+s)} Ei_1(t_j(a+s)) \qquad (26) $$

then,

$$ \mathcal{L}\{\eta(t)\} = \frac{C'_{0}}{s+a} + \sum_{j=1}^{N} a_j t_j e^{t_j(a+s)} Ei_1(t_j(a+s)) \qquad (27) $$

The corresponding Bode diagrams are given in Figure 8.

<figure><img src="image_10.png" /><figcaption>Figure 7. Fractional behaviour with a constant asymptote in low frequency.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8. Bode diagram of function (22).</figcaption></figure>

Kernel (23) is particularly interesting in a modelling context. Most of the phenomena modelled are of the low-pass filter type, which justifies the presence of the exponential term in (23) (the inverse Laplace of a low-pass filter is an exponential function). With this kernel,
the rational part in (23) is thus used to model the difference between the behaviour to be modelled (possibly fractional) and the low-pass behaviour provided by the exponential function. The next section shows the efficiency of kernels (4) and (23) to model fractional behaviour.

## 4. Kernels and Model Parameters Estimation

In this section, kernels (4) and (23) are used to model a fractional behaviour produced by a fractional transfer function supplied by a pseudo random binary sequence. We propose an algorithm to estimate the parameters used in the definition of the kernels with a nonlinear optimisation method. Other algorithms, especially for the computation of the kernel samples from a system input/output dataset using a least square method, will be presented in another paper.

### 4.1. Fractional Behaviour to Be Modelled

The data modelled in this part are produced by the fractional model:

$$H(s) = 5 \frac{\frac{s^{0.5}}{0.4} + 1}{\left(\frac{s^{0.5}}{0.01} + 1\right) \left(\frac{s^{0.5}}{2} + 1\right) \left(\frac{s^{0.5}}{100} + 1\right)} \quad (28)$$

The Bode diagram of transfer function (28) is shown in Figure 9 and highlights a fractional behaviour over a frequency band of five decades.

<figure><img src="image_12.png" /><figcaption>Figure 9. Frequency response of this model (28).</figcaption></figure>

To generate data that will be processed in the sequel, transfer function (28) is supplied by the pseudo random binary sequence of Figure 10a. Such a signal is currently used in a modelling context as its frequency content is almost constant over a wide frequency range. The time response, denoted $Y(t)$, is obtained through the approximation of the transfer function (28) using a nine poles and zeros distribution of each $s^{0.5}$ operator over the frequency range $[10^{-3} \text{ rd/s}, 10^3 \text{ rd/s}]$, leading to an approximation with a very large dimension: with 35 zeros and 36 poles. The transfer function output is represented by Figure 10b.
<figure><img src="image_13.png" /><figcaption>Figure 10. Pseudo random binary sequence (a) applied to model (28) and its time response (b).</figcaption></figure>

## 4.2. Nonlinear Optimization to Obtain Model Parameters

The input/output data set generated in Section 4.1 is then used to find the parameters of model (29) in which the kernel $\eta(t)$ is defined by (4) or (23). The goal is not to develop a robust identification methodology, but just to show that the kernels in Equations (4) and (23) are able to model a fractional behaviour.

The kernels

$$ \eta_{a_1}(t) = \frac{b_2 t^2 + b_1 t + 1}{c_2 t^2 + c_1 t + c_0}, \quad \eta_{a_2}(t) = \frac{b_2 t^2 + b_1 t + 1}{c_2 t^2 + c_1 t + c_0} e^{-at} \quad \text{and} \quad \eta_{a_3}(t) = \frac{b_3 t^3 + b_2 t^2 + b_1 t + 1}{c_3 t^3 + c_2 t^2 + c_1 t + c_0} \quad (29) $$

are thus considered. To find the set of five, six and seven parameters {$a, b_1, b_2, b_3, c_0, c_1, c_2, c_3$} depending on the kernel considered, a nonlinear Levenberg-Marquardt optimisation algorithm is used to minimize the criterion:

$$ \epsilon_{q_i} = \frac{1}{N_s} \sum_{k=1}^{N_s} (Y(kT_e) - Y_{estim_i}(kT_e))^2 \quad i \in \{1, 2, 3\}. \quad (30) $$

where $N_s$ is the number of samples and $Y_{estim_i}(kT_e)$, $i \in \{1, 2, 3\}$ is computed according to relation

$$ Y_{estim_i}(t) = \int_{0}^{t} \eta_{a_i}(t - \tau)u(\tau)d\tau \quad (31) $$

after integral discretisation in which the kernel involved is $\eta_{a_i}(t)$. The parameters and the criterion values obtained are given in Table 1. As the optimization domain defined by Equation (3) is not convex, there is a risk with this parameter estimation method that the minimization algorithm will become stuck on a local minimum.

<table><caption>Table 1. Parameters associated with kernels $\eta_{a_i}(t)$ and criterion values.</caption><thead><tr><th></th><th>a</th><th>b<sub>1</sub></th><th>b<sub>2</sub></th><th>b<sub>3</sub></th><th>c<sub>0</sub></th><th>c<sub>1</sub></th><th>c<sub>2</sub></th><th>c<sub>3</sub></th><th>ε<sub>qi</sub></th></tr></thead><tbody><tr><td>η<sub>a<sub>1</sub></sub>(t)</td><td></td><td>1.329 × 10<sup>2</sup></td><td>57.44</td><td></td><td>0.1121</td><td>1.604 × 10<sup>2</sup></td><td>4.827 × 10<sup>3</sup></td><td></td><td>3.952 × 10<sup>-7</sup></td></tr><tr><td>η<sub>a<sub>2</sub></sub>(t)</td><td>0.2018</td><td>2.188 × 10<sup>2</sup></td><td>2.046 × 10<sup>2</sup></td><td></td><td>0.1031</td><td>2.003 × 10<sup>2</sup></td><td>8.891 × 10<sup>3</sup></td><td></td><td>2.056 × 10<sup>-7</sup></td></tr><tr><td>η<sub>a<sub>3</sub></sub>(t)</td><td></td><td>3.683 × 10<sup>8</sup></td><td>8.747 × 10<sup>9</sup></td><td>2.247 × 10<sup>9</sup></td><td>0.0843</td><td>1.059 × 10<sup>8</sup></td><td>2.768 × 10<sup>10</sup></td><td>2.709 × 10<sup>11</sup></td><td>3.484 × 10<sup>-8</sup></td></tr></tbody></table>

The signals $Y_{estim_i}(t)$ and $Y(t)$ are compared in Figure 11, which also shows the modelling absolute error. This figure reveals the ability of kernels such as (4) or (23) to model fractional behaviours with a small number of parameters.
<figure><img src="image_14.png" /><figcaption>Figure 11. Comparison of $Y(t)$ and $Y_{estim_1}(t)$: (a) $Y_{estim_1}(t)$, (b) $Y_{estim_2}(t)$, (c) $Y_{estim_3}(t)$ and associated absolute error.</figcaption></figure>

5. Physical Interpretation and Realisation

Using Equation (16), an interpretation of the kernel in terms of delay distribution is possible. According to Appendix A, Equation (16) is indeed valid in the frequency domain, namely if $s = j\omega$. The following approximation is thus possible:

$$e^{j\omega b} \int_{\omega b}^{\infty} \frac{e^{-j\zeta}}{\zeta} d\zeta \approx e^{j\omega b} \sum_{k=0}^{\infty} \frac{e^{-j\omega b(1+k\Delta u)}}{1+k\Delta u} \quad (32)$$
This relation can be rewritten as:

$$e^{j\omega b} \int_{1}^{\infty} \frac{e^{-\zeta j\omega b}}{\zeta} d\zeta \approx \sum_{k=0}^{\infty} \frac{e^{-j\omega b k \Delta u}}{1 + k \Delta u}. \quad (33)$$

This result can be applied to Equation (15), thus making it possible to write:

$$\eta(j\omega) \approx \frac{C'_0}{j\omega} + \sum_{j=1}^{N} a_j t_j \sum_{k=0}^{\infty} \frac{e^{-j\omega t_j k \Delta u}}{1 + k \Delta u}. \quad (34)$$

This relation clearly highlights that the kernel $\eta(t)$ can be approximated by a series of delays, which can be compared to the analysis made in [34]. In [34], a probabilistic interpretation of fractional operators leads to their approximation in terms of delay distribution. It was shown that the approximation of fractional operators in terms of the distribution of exponential functions (or poles and zeros distribution) struggles to explain a large number of physical phenomena producing fractional behaviours. In these phenomena, such as diffusion, adsorption or aggregation, the fractional behaviour is produced by the stochastic behaviour of a large number of agents (atoms, molecules, people, etc.) that evolve in a constrained space to achieve a common goal (cross a material section in the case of diffusion, find a place on a surface in the case of adsorption, cling to a cluster in the case of aggregation, etc.). The delays characterise the time necessary for these agents to reach this objective and the coefficient which accompanies this delay characterises the density of the agents concerned.

This physical interpretation is now illustrated in the case of water diffusion. As shown by Figure 12, the medium is assumed to be constituted of channels of various lengths. Before the water hits the medium, the channels are assumed to be empty of fluid. When the fluid pressure $P(t)$ appears at the input side of the medium, it creates a flow inside the channels such that the total flow at the output side, denoted $Q(t)$, is the sum of the flow produced by each channel, denoted $Q_k(t)$:

$$Q(t) = \sum_{k} Q_{k}(t). \quad (35)$$

<figure><img src="image_15.png" /><figcaption>Figure 12. Interpretation in terms of time delay distribution of the fractional behaviour produced by diffusion phenomena.</figcaption></figure>

Each flow $Q_k(t)$ depends on the pressure $P(t)$ and on the inverse of the hydraulic resistance $R_k$ of each channel. Moreover, due to the difference in the length of the channels,
each flow $Q_k(t)$ reaches the output side with a time delay $T_k$ depending on the channel length, thus leading for each $Q_k(t)$ for the relation:

$$Q_k(t) = \frac{1}{R_k} P(t - T_k). \qquad (36)$$

and thus for the total flux:

$$Q(t) = \sum_k \frac{1}{R_k} P(t - T_k). \qquad (37)$$

If $Q_\infty$ denotes the steady state value of the flow $Q(t)$ when $P(t)$ is a unit step, then with such a modelling approach, it can be said that the coefficients $1/(Q_\infty R_k)$ can be viewed as the probability of having a delay in duration $T_k$ in the fractional behavior produced by the diffusion phenomena. More physically, $1/(Q_\infty R_k)$ can be connected to a distribution in channel length and to the probability of finding the corresponding channel in the studied system.

## 6. Conclusions

In this paper, fractional behaviours are modelled using convolution models involving a rational kernel. If the fractional behaviour is a pure power law function, a particular rational kernel made of recursively linked poles and zeros permits an accurate fitting on a time domain. This fitting solution is similar to the well-known one commonly used to approximate the behaviour of a fractional integrator in the frequency domain. But it is observed that the time domain solution here proposed produces less significant gain and phase oscillations for the same number of parameters (zeros and of poles). An identification method is also proposed in this paper to find the rational kernel parameters in the case of any fractional behaviour (not a pure power law function), making it possible to work directly in the time domain. It is applied on data generated by a fractional transfer function supplied by a pseudo random binary sequence. The interest in the proposed rational kernels is highlighted by this application. The rational kernels that we have used also permit a physical interpretation of the modelled phenomena. Using Laplace transform, it is shown that the rational kernels studied in the paper can be realised by a series of delays. This realisation is close to the probabilistic interpretation of the fractional behaviours made in [3,34], with a physical meaning in relation to phenomena such as diffusion, adsorption and aggregation, which are known to produce fractional behaviours. This work is an additional illustration that fractional behaviours can be modelled using tools other than those involving fractional calculus.

The authors also have to now consider additional real-life fractional behaviour in order to clearly identify the limits of the use of non-singular rational kernels in modelling situations and to consider possible significant changes over time in the modelled system behaviour or possible nonlinearities. The authors also intend to develop other parameter estimation algorithms and thus discretisation strategies [35]. We also want to develop tools to analyse the properties of the convolution models (such as the one given by Equation (31)) involving rational non-singular kernels in terms of stability, controllability and observability.

**Author Contributions:** Conceptualization, J.S. and C.F.; methodology, J.S. and C.F.; writing—original draft preparation, J.S. and C.F.; writing—review and editing, J.S. and C.F.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data are available on demand.

**Conflicts of Interest:** The authors declare no conflict of interest.
## Appendix A. Appendix Evaluation of the Exponential Integral on the Imaginary Axis

According to [32], the exponential integral function is defined by [33]:

$$Ei_1(x) = \int_x^\infty \frac{e^{-\zeta}}{\zeta} d\zeta \quad \text{with} \quad |\arg x| < \pi \quad (A1)$$

where $x = z + j\omega$. The integral converges if the upper limit is $\infty e^{j\alpha}$ [36], and is independent of $\alpha$ as long as $-\frac{\pi}{2} \le \alpha \le \frac{\pi}{2}$. Let $\zeta = j\tau$, and then Equation (A1) becomes

$$Ei_1(z + j\omega) = \int_{\omega-jz}^{\infty} \frac{e^{-j\tau}}{\tau} d\tau, \quad (A2)$$

and thus

$$Ei_1(z + j\omega) = \int_{\omega-jz}^{\infty} \frac{\cos(\tau)}{\tau} d\tau - j \int_{\omega-jz}^{\infty} \frac{\sin(\tau)}{\tau} d\tau. \quad (A3)$$

Introducing sine and cosine integral function with $|\arg x| < \pi$, respectively, is defined by [33]:

$$Ci(x) = \int_x^\infty \frac{\cos(\tau)}{\tau} d\tau \quad \text{and} \quad Si(x) = \int_0^x \frac{\sin(\tau)}{\tau} d\tau \quad \text{with} \quad -\int_x^\infty \frac{\sin(\tau)}{\tau} d\tau = Si(x) - \frac{\pi}{2} \quad (A4)$$

Equation (A2) becomes for $z = 0$

$$Ei_1(j\omega) = -Ci(\omega) + j\left(Si(\omega) - \frac{\pi}{2}\right). \quad (A5)$$

This last relation permits us to evaluate the frequency resulting from of $Ei_1(s)$, with s being the Laplace variable.

For $z = 0$, Equation (A3) also becomes the following:

$$Ei_1(j\omega) = \int_{\omega}^{\infty} \frac{\cos(\tau)}{\tau} d\tau - j \int_{\omega}^{\infty} \frac{\sin(\tau)}{\tau} d\tau = \int_{\omega}^{\infty} \frac{e^{-j\tau}}{\tau} d\tau \quad (A6)$$

## References

1. Podlubny, I. *Fractional Differential Equations: Theoretical Developments and Applications*. In *Physics and Engineering Mathematics in Sciences and Engineering*; Academic Press: Cambridge, MA, USA, 1999.

2. Sabatier, J.; Agrawal, O.P.; Machado, J.A.T. *Advances in Fractional Calculus*; Springer: Dordrecht, The Netherlands, 2007.

3. Sabatier, J.; Farges, C.; Tartaglione, V. *Fractional Behaviours Modelling: Analysis and Application of Several Unusual Tools*; Intelligent Systems, Control and Automation: Science and Engineering Series; Springer: Cham, Switzerland, 2022; Volume 101.

4. Samko, S.G.; Kilbas, A.A.; Marichev, O.I. *Fractional Integrals and Derivatives: Theory and Applications*; Gordon and Breach Science Publishers: London, UK, 1993.

5. Monje, C.A.; Chen, Y.; Vinagre, B.M.; Xue, D.; Feliu, V. *Fractional-Order Systems and Controls: Fundamentals and Applications*; Advances in Industrial Control Series; Springer: London, UK, 2010.

6. Baleanu, D.; Balas, V.E.; Agarwal, P. *Fractional Order Systems and Applications in Engineering*; Elsevier: Amsterdam, The Netherlands, 2022.

7. Petras, I. *Fractional Order Systems*; MDPI Books; MDPI: Basel, Switzerland, 2019.

8. Caputo, M.; Fabrizio, M. A new definition of fractional derivative without singular kernel. *Prog. Fract. Differ. Appl.* **2015**, *1*, 73–85.

9. Dokoumetzidis, A.; Magin, R.; Macheras, P. A commentary on fractionalization of multi-compartmental models. *Pharm. Pharm.* **2010**, *37*, 203–207. [CrossRef] [PubMed]

10. Sabatier, J.; Farges, C.; Trigeassou, J.C. Fractional systems state space description: Some wrong ideas and proposed solutions. *J. Vib. Control* **2014**, *20*, 1076–1084. [CrossRef]

11. Balint, A.M.; Balin, S. Mathematical description of the groundwater flow and that of the impurity spread which use temporal Caputo or Riemann-Liouville fractional partial derivatives is non-objective. *Fractal Fract.* **2020**, *4*, 36. [CrossRef]

12. Pantokratoras, A. Comment on the paper “Fractional order model of thermo-solutal and magnetic nanoparticles transport for drug delivery applications, Subrata Maiti, Sachin Shaw, G.C. Shit, [Colloids Surf. B Biointerfaces, 203(2021) 111754]”. *Colloids Surf. B Biointerfaces* **2023**, *222*, 113074. [CrossRef]
13. Pantokratoras, A. Discussion on the paper "A Numerical Scheme for Fractional Mixed Convection Flow Over Flat and Oscillatory Plates, Yasir Nawaz, Muhammad Shoaib Arif, Kamaleldin Abodayeh, Journal of Computational and Nonlinear Dynamics, July 2022, Vol. 17, 071008". J. Comput. Nonlinear Dyn. **2023**, *18*, 115501. [CrossRef]

14. Atangana, A.; Baleanu, D. New fractional derivatives with nonlocal and non-singular kernel: Theory and application to heat transfer model. *Therm. Sci.* **2016**, *20*, 763–769. [CrossRef]

15. Hattaf, K. A new generalized definition of fractional derivative with non-singular kernel. *Computation* **2020**, *8*, 49. [CrossRef]

16. Hristov, J. Response functions in linear viscoelastic constitutive equations and related fractional operators. *Math. Model. Nat. Phenom.* **2019**, *14*, 305. [CrossRef]

17. Fernandez, A.; Al-Refai, M. A Rigorous Analysis of Integro-Differential Operators with Non-Singular Kernels. *Fractal Fract.* **2023**, *7*, 213. [CrossRef]

18. Sun, H.; Wang, Y.; Yu, L.; Yu, X. A discussion on nonlocality: From fractional derivative model to peridynamic model. *Commun. Nonlinear Sci. Numer. Simul.* **2022**, *114*, 106604. [CrossRef]

19. Sabatier, J. Modelling Fractional Behaviours Without Fractional Models. *Front. Control. Eng.* **2021**, *2*, 716110. [CrossRef]

20. Newman, D.J. Rational approximation to $|x|$. *Mich. Math. J.* **1964**, *11*, 11–14. [CrossRef]

21. Manabe, S. The non-integer Integral and its Application to control systems. *ETJ Jpn.* **1961**, *6*, 83–87.

22. Carlson, G.E.; Halijak, C.A. Simulation of the Fractional Derivative Operator and the Fractional Integral Operator. Available online: http://krex.k-state.edu/dspace/handle/2097/16007 (accessed on 18 September 2023).

23. Ichise, M.; Nagayanagi, Y.; Kojima, T. An analog simulation of non-integer order transfer functions for analysis of electrode processes. *J. Electroanal. Chem. Interfacial Electrochem.* **1971**, *33*, 253–265. [CrossRef]

24. Oustaloup, A. *Systèmes Asservis Linéaires D'ordre Fractionnaire*; Masson: Paris, France, 1983.

25. Raynaud, H.F.; Zergaïnoh, A. State-space representation for fractional order controllers. *Automatica* **2000**, *36*, 1017–1021. [CrossRef]

26. Charef, A. Analogue realisation of fractional-order integrator, differentiator and fractional PI$\lambda$D$\mu$ controller. *IEE Proc. Control Theory Appl.* **2006**, *153*, 714–720. [CrossRef]

27. Sabatier, J. Solutions to the Sub-Optimality and Stability Issues of Recursive Pole and Zero Distribution Algorithms for the Approximation of Fractional Order Models. *Algorithms* **2018**, *11*, 103. [CrossRef]

28. Kochubei, A.N. General fractional calculus, evolution equations, and renewal processes. *Integral Equ. Oper. Theory* **2011**, *71*, 583–600. [CrossRef]

29. Luchko, Y.; Yamamoto, M. The general fractional derivative and related fractional differential equations. *Mathematics* **2020**, *8*, 2115. [CrossRef]

30. Luchko, Y. General fractional integrals and derivatives and their applications. *Phys. D Nonlinear Phenom.* **2023**, *455*, 133906. [CrossRef]

31. Tarasov, V.E. Multi-kernel general fractional calculus of arbitrary order. *Mathematics* **2023**, *11*, 1726. [CrossRef]

32. Tseng, P.-H.; Lee, T.-C. Numerical evaluation of exponential integral: Theis well function approximation. *J. Hydrol.* **1998**, *205*, 38–51. [CrossRef]

33. Abramowitz, M.; Stegun, I.A. *Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables*; National Bureau of Standards, Applied Mathematics Series, No. 55; U.S. Government Printing Office: Washington, DC, USA, 1964; first printing.

34. Sabatier, J. Probabilistic Interpretations of Fractional Operators and Fractional Behaviours: Extensions, Applications and Tribute to Prof. José Tenreiro Machado's Ideas. *Mathematics* **2022**, *10*, 4184. [CrossRef]

35. Ortigueira, M.D.; Machado, J.T. The 21st century systems: An updated vision of discrete-time fractional models. *IEEE Circuits Syst. Mag.* **2022**, *22*, 36–56. [CrossRef]

36. Corrington, M.S. Applications of the Complex Exponential Integral. *Math. Comput.* **1961**, *15*, 1–6. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.