# Contribution to the initialization of linear non-commensurate fractional-
order systems for the joint estimation of parameters and fractional
differentiation orders

Mohamed A. Bahloul1,2, Zehor Belkhatir3 and Taous-Meriem Laleg Kirati2,4
1 College of Engineering, Electrical Engineering Department at Alfaisal
University,
Riyadh 11533, Saudi Arabia. E-mail<EMAIL_ADDRESS>
2 Electrical and Computer Engineering Department, KAUST, Saudi Arabia.
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
3 School of Engineering and Sustainable Development, De Montfort University,
United Kingdom, E-mail<EMAIL_ADDRESS>
4 National Institute for Research in Digital Science and Technology, Paris-
Saclay, France

###### Abstract

It has been recognized that using time-varying initialization functions to
solve the initial value problem of fractional-order systems (FOS) is both
complex and essential in defining the dynamical behavior of the states of
FOSs. In this paper, we investigate the use of the initialization functions
for the purpose of estimating unknown parameters of linear non-commensurate
FOSs. In particular, we propose a novel "pre-initial" process that describes
the dynamic characteristic of FOSs before the initial state and consists of
designing an appropriate time-varying initialization function that ensures
accurate convergence of the estimates of the unknown parameters. To do so, we
propose an estimation technique that consists of two steps: (i) to design of
practical initialization function that is output-dependent and which is
employed; (ii) to solve the joint estimation problem of both parameters and
fractional differentiation orders (FDOs). A convergence proof has been
presented. The performance of the proposed method is illustrated through
different numerical examples. Potential applications of the algorithm to joint
estimation of parameters and FDOs of the fractional arterial Windkessel and
neurovascular models are also presented using both synthetic and real data.
The added value of the proposed "pre-initial" process to solve the studied
estimation problem is shown through different simulation tests that
investigate the sensitivity of estimation results using different time-varying
initialization functions.

##  1 Introduction

Over the past 325 years, fractional calculus (FC) has attracted the attention
of mathematicians and engineers working in various fields of science and
engineering [1]. FC started to be used as a powerful tool to describe and
explore complicated dynamical systems of real-world applications, e.g.,
biology [2], control [3], economic [4]. Indeed FC can describe systems with
high-order dynamics and complex nonlinear aspects employing fewer coefficients
than the integer-order calculus. The arbitrary order of the derivatives
provides an extra degree of freedom to concisely and precisely entail the
hidden dynamics having a different origin of memory effect. A distinctive
feature of fractional-order calculus is that contrary to integer-order
calculus, the fractional-order one accounts for the system’s memory. In fact,
the fractional-order derivative depends not only on the local conditions of
the evaluated time but also on the entire history of the function. This
peculiarity is usually valuable when the studied system holds a long-term
"memory," and any assessed point depends on the past values of the function
[5]. Accordingly, one essentially-crucial prerequisite to ensuring the
fractional-order theory functioning is the adequate initialization and the
proper incorporation of the history of the system [6].

Generally, the initialization of fractional-order systems is considered
complicated and challenging. This dilemma surfaced notably in the case of the
non-zero initial value. In fact, it is challenging to design appropriate
initialization functions that guarantee the acquiring of exact states of the
system while accounting for its history, [7]. The pre-initialization process
and the initialization functions of non-zero initial value-based fractional-
order systems remain an open and controversial question. Various effective
methods have been developed to analyze the characteristic of fractional order
derivatives where initial values are insufficient [8, 9, 10]. Lorenzo et al.,
[11] was the first to demonstrate that time-varying functions are more
suitable for the initialization of fractional-order systems rather than
constant ones. The time-varying initialization has a profound effect on the
standard definitions of fractional-order derivative and integral. Different
time-varying initialization functions may lead to the same initial value from
where the fractional-order operator starts; however, as the system’s dynamic
is related to the pre-initial process, different initialization functions
result in different responses. This fact is known as the aberration phenomenon
[12]. Accordingly, designing the proper initialization function to acquire the
desired convergence of the fractional differential system is deemed very
problematic.

For this reason, in this report, we introduce a novel pre-initialization
process that warrants a fast and precise convergence of the joint estimation
of the parameters and differentiation orders of the fractional differential
system (FOS). The key hypothesis in the design is to consider an output-
dependent initialization function when estimating the unknowns. This will
reduce the infinite-dimensional space of initialization functions into a
finite parametric space where the remaining degree of freedom is about the
length of the history function to consider. In addition, we proposed a two-
stage algorithm: while the first stage solves a system of linear equations to
estimate the parameters of the FOS, the second stage uses the iterative first-
order Newton’s method [13, 14, 15, 16] to estimate the fractional
differentiation orders. Our contributions are as follows: (i) we consider a
time-varying function to initialize the FOS; (ii) we design a pre-
initialization process based on the output signal; and (iii) we solve a simple
linear equation system to estimate the parameters, and iteratively we apply
Newton’s method to estimate the fractional differentiation orders.

The performance of the proposed method is illustrated through different
numerical examples. Additionally, potential applications of the algorithm are
presented, which consists of estimating parameters and fractional
differentiation orders of a fractional-order arterial Windkessel and
neurovascular models. To the best to the author’s knowledge, this is the first
study that accounts for the initialization function as part of the parameters
estimation problem for fractional differential systems.

##  2 Notations

For a clear perception of the report, in this section, we present the adopted
notation and the basic definitions of the fractional-order integral and
derivative. Through the following we consider a smooth function $f(t)$ such
that $f(t)$ is zero for $t\leq t_{abs}$ and $f(t)$ is $f_{in}(t)$ for
$t_{abs}\leq t\leq t_{in}$.

  * •

${}_{t_{in}}D_{t}^{\alpha}f(t)$ denotes the ’initialized’ $\alpha^{th}$ order
differ-integration of $f(t)$ from start point $t_{0}$ to $t$.

  * •

${}_{t_{in}}d_{t}^{\alpha}f(t)$ represents the ’non-initialized’ generalized
(or fractional) $\alpha^{th}$ order differ-integration of $f(t)$. It is
equivalent to shifting the origin of function $f(t)$ at the start of the point
from where differ-integration starts.

$\dfrac{d^{\alpha}f(t)}{[d(t-{t_{in}})^{\alpha}]}\equiv\ \
_{{t_{in}}}d_{t}^{\alpha}f(t).$ (1)

  * •

$\Psi$ denotes the time-varying initialization function or history function.
It is also known as the complementary function [17]. The expression between
initialized differ-integral and non-initialized ones is:

${}_{t_{in}}D_{t}^{\alpha}f(t)=_{t_{in}}d_{t}^{\alpha}f(t)+\Psi(f,\alpha,t_{abs},t_{in},t)$
(2)

###  2.1 Definitions

There are many kinds of definitions for fractional-order differ-integration.
Here we present the more commonly known ones, namely the Riemann-Liouville,
Caputo and Grunwald-Letnikov definitions. A detailed note on the different
definitions might be found in this reference, [11, 17].

######  Definition 1.

The Riemann-Liouville (R-L) definition of fractional-order integration is
drawn as:

$_{t_{in}}^{RL}d_{t}^{-\alpha}f\left(t\right)=\frac{1}{\Gamma\left(\alpha\right)}\int_{t_{in}}^{t}\frac{f\left(\tau\right)}{\left(t-\tau\right)^{\alpha}}d\tau,$
(3)

here $t_{abs}$ is noted as the terminal point as well. $0<\alpha<1$ and
$\Gamma(x)$ corresponds to so-called Gamma-function, written as
$\Gamma(x)=\int_{0}^{\infty}e^{-u}u^{x-1}du.$

######  Definition 2.

The R-L definition of fractional-order derivative is based on the above
definition (3) and the standard integer-order derivative:

$_{t_{in}}^{RL}d_{t}^{\alpha}f\left(t\right)=\frac{\mathrm{d}}{\mathrm{d}t}\left[{}_{t_{in}}d_{t}^{-(1-\alpha)}f\left(t\right)\right].$
(4)

######  Definition 3.

The Caputo definition of fractional-order differentiation takes the integer-
order differentiation of the function first and then takes a fractional-order
integration:

$_{t_{in}}^{C}d_{t}^{\alpha}f\left(t\right)=\dfrac{1}{\Gamma\left(n-\alpha\right)}\int_{{t_{in}}}^{t}\left(1-\tau\right)^{-\alpha-1+n}\dfrac{\mathrm{d^{n}}}{\mathrm{d}t^{n}}f\left(\tau\right)d\tau.$
(5)

######  Definition 4.

The Grunwald-Letnikov based Fractional-order differ-integrationG-L defines the
fractional integration and differentiation in a unified way:

$_{t_{in}}^{GL}d_{t}^{\alpha}f\left(t\right)=\lim_{h\to
0}h^{-\alpha}{\sum}_{j=0}^{{\left\lfloor\dfrac{t-t_{in}}{T_{s}}\right\rfloor}}(-1)^{j}\binom{\alpha}{j}f(t-j\,\,h),$
(6)

where $h$ is the sampling time,
${\binom{\alpha}{j}}=\dfrac{\Gamma(\alpha+1)}{\Gamma(j+1)\Gamma(\alpha-j+1)}$
and $\left\lfloor.\right\rfloor$ describes the floor function and denotes the
biggest integer smaller or equal to the argument.

######  Remark 1.

Suppose $f\in C^{\left\lfloor\alpha\right\rfloor+1}$, the G-L derivative is
equivalent to the R-L derivative (see [18]).

###  2.2 Aberration phenomenon

This section focuses on the proposition of the so-called aberration phenomena
that represent the effect of the pre-initial processes on a fractional-order
system’s response. To do so, we simulate a simple system that model the
internal force of an axially loaded viscoelastic bar, [19]. In the following
illustration example, $u(t)$ is the input denoting the longitudinal force, and
$y(t)$ is the output representing the elongation. They satisfy the following
fractional-order relationship:

${}_{t_{0}}D_{t}^{0.5}y(t)=u(t)$ (7)

Assume the force $u(t)$, the input of the system, is zero before the absolute
zero instant of the system, and the initial instant is taken as $t_{in}=5$. We
excite the system by three different inputs $U_{0}$, $U_{1}$, and $U_{2}$ that
are specified as follows:

$\begin{matrix}U_{0}\left\\{\begin{matrix}0,\ \ t\leqslant 0\\\ 0,\ \ 0<t<5\\\
1,t\geqslant 5\end{matrix}\right.,&U_{1}\left\\{\begin{matrix}0,\ \ t\leqslant
0\\\ t,\ \ 0<t<5\\\ 1,t\geqslant
5\end{matrix}\right.,&U_{2}\left\\{\begin{matrix}0,\ \ t\leqslant 0\\\
\dfrac{1}{4}t^{4},\ \ 0<t<5\\\ 1,t\geqslant 5\end{matrix}\right.\end{matrix}$
(8)

Figure 1:  The aberration phenomenon.

The corresponding outputs ($Y_{0}$, $Y_{1}$, and $Y_{2}$) along with the
initialization function, as developed in the previous section, are shown in
Figure 1. It is clear from this graph that although $Y_{1}$ and $Y_{2}$ have
the same values at the initial time $t_{in}=5$ and they are excited with the
same input afterward, they don’t have the same response later. In addition,
from the output $Y_{0}$ that corresponds to a non-initialized response as
$u(t)$ is equal to zero over the pre-initialized process $t\in[0\ \ 5]$. From
all these outputs, we can perceive the effect of the pre-initial process.
Accordingly, in the simulation of the fractional-order system, the ignorance
of infinite states at $t_{in}$ or the improper initialization of the states
before the initial instant might cause an aberration of the output. On the one
hand, the noticed aberration phenomenon reflects the complexity and the
difficulty of the initial value problem. On the other hand, it illustrates the
great challenge of analyzing and controlling fractional-order systems.
Commonly, to overcome the negative consequences of the aberration, former
studies are used to transform the problems of the non-initial value into zero
initial value ones. Hence, it is essential to seek the inner essence of this
aspect. Furthermore, it is primordial to provide adapted procedures to select
the exact time-varying history function that leads to the correct response and
accurate estimation of the system’s states.

###  2.3 Time-varying function based-initialization of FOS

The main idea of the history function based-initialization is that the
initialization time-varying function $\Psi(t)$ should bring out the past
history. Otherwise stated, the history function entails the effect of
fractionally integrating the function from its birth [17, 15].
Assume that $f(t)$ was born at $t=t_{abs}$, that is $f(t)=0$ for all time less
than equal to $t_{abs}$. Then the time period between $t_{abs}$ and $t_{in}$
may be considered as history if the fractional-order integration starts at
$t=t_{in}$ and $f(t)=f_{in}(t)$ $\forall\,\,\,t\in[t_{abs}\,\,\,t_{in}]$. The
fundamental idea is that the fractional-order integral
(${}_{t_{in}}d_{t}^{-\alpha}f(t)$), should be properly initialized so that it
should function as continuation of integral starting at $t=t_{abs}$. Hence an
initialization function, $\Psi$ must be added to
(${}_{t_{in}}d_{t}^{-\alpha}f(t)$), so that the fractional-order integral
starting at $t=t_{in}$ should be identical to the result starting at
$t=t_{abs}$ for $t\geq t_{in}$. The history function $\Psi$ has the effect of
allowing the function $f(t)$ and its derivatives to start at a value other
than $t_{in}$, namely ${}_{t_{abs}}D_{t_{in}}^{-\alpha}f(t)\mid_{t=t_{in}}$ ,
and continues to contribute to differ-integral response after $t=t_{in}$. That
is, a function of time is added to the uninitialized integral, (not just a
constant at $t=t_{in}$). The above argument can be formulated as follows: for
$t>t_{in}$

$\Psi(f_{in},\alpha,t_{abs},t_{in},t)=_{t_{abs}}d_{t}^{-\alpha}f(t)-_{t_{in}}d_{t}^{-\alpha}f(t).$
(9)

Generally, we differentiate two types of initialization:

  * •

Side-initialization: Fully arbitrary initialization may be applied to the
differ-integral operator at time $t=t_{in}$ .

  * •

Terminal initialization: It is assumed that the differ-integral operator can
be initialized (charged) by effectively differ-integrating prior to the start
time, $t=t_{in}$.

Figure 2 demonstrates the concept of initialization as a block diagram of a
signal flow graph.

Figure 2:  Signal flow graph for demonstrating initialization of fractional
integration [17]. (a) illustrates the Side-initialization and (b) Terminal
initialization

####  2.3.1 History-function based-initialization of fractional-order
integration

Assume, the fractional-order integration of a function $f(t)$ starts a $t=c$
and takes place for all $t>t_{in}\geqslant t_{abs}$. The initialization
process takes place during the period $t\in[t_{abs}\ \ t_{in}]$. Generally the
Riemann-Liouville (R-L) definition is admitted in the case where the
differintegrand $f(t)=0$ for all $t\leqslant t_{abs}$. The fractional-order
integration based R-L definition of $f(t)$ is:

${}_{t_{abs}}D_{t}^{-\alpha}f(t)=_{t_{abs}}d_{t}^{-\alpha}f(t)=\dfrac{1}{\Gamma(\alpha)}\int_{t_{abs}}^{t}(t-\tau)^{\alpha-1}f(\tau)d\tau,$
(10)

where ($\alpha\geqslant 0$ and $t>t_{abs}$ subject to $f(t)=0$ for all
$t\leqslant t_{abs}$. $\Psi(f,-\alpha,t_{abs},t_{abs},t)=0$
The following definition of fractional-order integration will apply generally
(at any $t>c$ ):

${}_{c}D_{t}^{-\alpha}f(t)=\dfrac{1}{\Gamma(\alpha)}\int_{c}^{t}(t-\tau)^{\alpha-1}f(\tau)d\tau+\Psi(f,-\alpha,a,c,t)$
(11)

where ($\alpha\geqslant 0$, $t>t_{abs}$, $t_{in}>t_{abs}$ and $f(t)=0$ for all
$t\leqslant t_{abs}$.
The function $\Psi(f,-\alpha,t_{abs},t_{in},t)$ is called the initialization
function and will be chosen such that

${}_{t_{abs}}D_{t}^{-\alpha}f(t)=_{t_{in}}D_{t}^{-\alpha}f(t),\ \ t>t_{in}$
(12)

The above condition gives:

$\dfrac{1}{\Gamma(\alpha)}\int_{t_{abs}}^{t}(t-\tau)^{\alpha-1}f(\tau)d\tau=\dfrac{1}{\Gamma(\alpha)}\int_{t_{in}}^{t}(t-\tau)^{\alpha-1}f(\tau)d\tau+\Psi(f,-\alpha,t_{abs},t_{in},t)$
(13)

Hence

$\begin{matrix}\Psi(f,-\alpha,t_{abs},t_{in},t)=\dfrac{1}{\Gamma(\alpha)}\int_{t_{abs}}^{t}\left(t-\tau)^{\alpha-1}f(\tau\right)d\tau-\dfrac{1}{\Gamma(\alpha)}\int_{t_{in}}^{t}(t-\tau)^{\alpha-1}f(\tau)d\tau\end{matrix}$
(14)

Since
$\int_{t_{abs}}^{t}g(\tau)d\tau=\int_{t_{abs}}^{t_{in}}g(\tau)d\tau+\int_{t_{in}}^{t}g(\tau)d\tau$
we get:

$\Psi(f,-\alpha,t_{abs},t_{in},t)=_{t_{abs}}D_{t_{in}}^{-\alpha}=\dfrac{1}{\Gamma(\alpha)}\int_{t_{abs}}^{t_{in}}\left(t-\tau)^{\alpha-1}f(\tau\right)d\tau$
(15)

This expression for $\Psi(t)$ gives ’terminal initialization’, and also brings
out in the definition of fractional integral the effect of the past ’history’,
namely the effect of fractionally integrating the $f(t)$ from $t_{abs}$ to
$t_{in}$. This effect is also called terminal charging.

####  2.3.2 History-function based-initialization of fractional-order
derivative

Unlike the integer-order derivative that is considered a local quantity, the
fractional-order derivative (FD) is a non-local operator and has a history.
Furthermore, as FD contains fractional-order integration, its evaluation
requires the initialization process. By considering the integer-order
derivative as a particular case of the fractional-order derivative, a
generalization of the integer-order derivative concept also calls for
initialization. Accordingly, a generalized integer-order differentiation with
initialization can be defined as:

${}_{t_{in}}D_{t}^{m}f(t)=\dfrac{d^{m}}{dt^{m}}f(t)+\Psi(f,m,t_{abs},t_{in},t),t>t_{in}$
(16)

where $m$ is an integer corresponding to the differentiation order and
$\Psi(f,m,t_{abs},t_{in},t)$ is the initialization function. Assume a
fractional-order differentiation; $\alpha\geqslant 0$. A function $f(t)$ is
born at $t=t_{abs}$ and before that the value is zero ($f(t)=0$ for
$t\geqslant t_{abs}$). $\alpha=m-\beta$ equivalently $m$ is the integer just
greater than the fractional order $\alpha$, by amount $\beta$. The
differentiation starts at $t>t_{in}$.
A non-initialized fractional-order derivative can be expressed as:

${}_{t_{abs}}d_{t}^{\alpha}f(t)=_{t_{abs}}D_{t}^{m}\
{}_{t_{abs}}D_{t}^{-\beta}f(t)$ (17)

Similar to the fractional-order integration in this case the initialization
function $\Psi(f,-\beta,t_{abs},t_{abs},t)=0$.
If we consider the $h(t)=_{t_{abs}}D_{t}^{-\beta}f(t)$ i.e. fractional
integral of function starting at $t_{abs}$ with initialized term
$\Psi(h,m,t_{abs},t_{abs},t)=0$; the initialized fractional derivative looks
and defined as for $\alpha\geqslant 0$, $t>t_{in}\geqslant t_{abs}$.

${}_{t_{in}}D_{t}^{\alpha}f(t)=_{t_{in}}D_{t}^{m}f(t)\
_{t_{in}}D_{t}^{-\beta}f(t)$ (18)

Accordingly, the history-function based-initialization of the fractional-order
derivative is generally similar to that obtained with initialization of the
fractional-order integral. In fact, it requires the following as the
fractional-order integral:

Figure 3:  Initialization of fractional derivative.

${}_{t_{abs}}D_{t}^{\alpha}f(t)=_{t_{in}}D_{t}^{\alpha}f(t)$ (19)

For all $t>t_{in}\geqslant t_{abs}$. Specifically this requires compatibility
of the derivatives starting at $t=t_{abs}$ and $t=t_{in}$ for all $t>t_{in}$.
Therefore it follows that:

${}_{t_{in}}D_{t}^{m}\ {}_{t_{in}}D_{t}^{-\beta}f(t)=_{t_{abs}}D_{t}^{m}\
{}_{t_{abs}}D_{t}^{-\beta}f(t)$ (20)

Expanding the fractional-order integrals with the initialization function
introduced in the previous section we obtain:

$\begin{matrix}{}_{t_{in}}D_{t}^{m}\left(\dfrac{1}{\Gamma(\beta)}\int_{t_{in}}^{t}(t-\tau)^{\beta-1}f(\tau)d\tau+\Psi(f,-\beta,t_{abs},t_{in},t)\right)=\\\
{}_{t_{abs}}D_{t}^{m}\left(\dfrac{1}{\Gamma(\beta)}\int_{t_{abs}}^{t}(t-\tau)^{\beta-1}f(\tau)d\tau+\Psi(f,-\beta,t_{abs},t_{abs},t)\right)\end{matrix}$
(21)

For, $t>t_{in}$ and $\Psi(f,-\beta,t_{abs},t_{abs},t)=0$. Using the definition
of generalized integer order derivative as defined above we get:

$\begin{matrix}\dfrac{d^{m}}{dt^{m}}\left\\{\dfrac{1}{\Gamma(\beta)}\int_{t_{in}}^{t}(t-\tau)^{\beta-1}f(\tau)d\tau+\Psi(f,-\beta,t_{abs},t_{in},t)\right\\}+\Psi(h_{1},m,t_{abs},t_{in},t)=\\\
\dfrac{d^{m}}{dt^{m}}\dfrac{1}{\Gamma(\beta)}\int_{t_{abs}}^{t}(t-\tau)^{\beta-1}f(\tau)d\tau+\Psi(h_{2},m,t_{abs},t_{abs},t)\end{matrix}$
(22)

where $h_{1}=_{t_{in}}D_{t}^{-\beta}f(t)$ and
$h_{2}=_{t_{abs}}D_{t}^{-\beta}f(t)$. The integer order derivative is
initialized at $t=a$ thus $\Psi(h_{2},m,t_{abs},t_{abs},t)=0$. After
rearranging the integrals we get:

$\Psi(h_{1},m,t_{abs},t_{in},t)=\dfrac{d^{m}}{dt^{m}}\left(\dfrac{1}{\Gamma(\beta)}\int_{t_{abs}}^{t_{in}}(t-\tau)^{\beta-1}f(\tau)d\tau-\Psi(f,-\beta,t_{abs},t_{in},t)\right)$
(23)

The above expression is the requirement for the initialization of the
derivative in general. Under the condition of terminal charging of the
fractional-order integral, the initialization function of fractional
integration as defined and derived earlier is:

$\Psi(f,-\beta,t_{abs},t_{in},t)=\dfrac{1}{\Gamma(\beta)}\int_{t_{abs}}^{t_{in}}(t-\tau)^{\beta-1}f(\tau)d\tau$
(24)

Hence the $\Psi(h_{1},m,t_{abs},t_{in},t)=0$. Figure 3 illustrates the
initialization concept for fractional derivative.

####  2.3.3 Grunwald-Letnikov based fractional-order differ-integrals
initialization concept

######  Definition 5.

The formulations of the time-varying history function using the most known
fractional-order differ-integration definitions namely RL, C and GL can be
observed in [11]. Here we present only the G-L based-one as our proposed
method is based on it.

$\begin{matrix}\Psi_{GL}(f_{in},\alpha,t_{in},t_{abs},t)=\underset{{h\rightarrow
0}}{\lim}\left\\{\dfrac{h^{-\alpha}}{\Gamma(-\alpha)}{\sum}_{j=0}^{N_{3}-1}\dfrac{\Gamma\left(N_{1}-1-\alpha-j\right)}{\Gamma\left(N_{1}-j\right)}f[t-(N_{1}-1-j)h]\\!\\!\right\\},\end{matrix}$
(25)

where ${h={(t-t_{abs})}\setminus{N_{1}}}$ is the time step, $t>t_{in}>t_{abs}$
and $N_{1}$ and $N_{3}$ are integers such that:

${N_{3}=\left(\dfrac{t_{in}-t_{abs}}{t-t_{abs}}\right)N_{1}}$ (26)

##  3 Problem Formulation

Let us consider the following FOS:

$y(t)+{\sum}^{N}_{i=1}a_{i}\,\,\,\,{}_{t_{in}}D^{\alpha_{i}}_{t}y(t)=b\,\,u(t),\,\,\,\,\,\,\,\,\,\,\,\,\,\,t\in[t_{in}\,\,\,\,\,\,T]$
(27)

where $y(t):[t_{in}\,\,\,T]\rightarrow\mathbb{R}$ is the output,
$u(t):\,\,\,[t_{in}\,\,\,T]\rightarrow\mathbb{R}$ is the input. $(a_{1}\
a_{2}\ \dots a_{N})$ and $b$ are real parameters.
$\alpha_{i}\in\wp=(n_{i-1},n_{i})$
111$\wp=\prod^{N}_{i=1}\wp_{i}=\\{(a_{1},a_{2}\dots,a_{N})|\forall
i=1,2,\dots,N,a_{i}\in\wp_{i}\\}$ is the generalized Cartesian product of $N$
sets $\wp_{1},\wp_{2}\dots\wp_{N}$. , with $n_{i}\in\mathbb{N}^{*}$ and
$i=1,2,\dots,N$ are the fractional differentiation orders. They are assumed to
be as follows:
${0\leq\\!\\!\alpha_{1}\leq\\!\\!\alpha_{2}\leq\dots\leq\\!\\!\alpha_{n}}$,
i.e., ${n_{i}<n_{i+1}}$ for ${i=1,2,\dots,N-1}$. We substitute the initialized
fractional derivative (2) in (27) which leads the following equation:

$y(t)+{\sum}^{N}_{i=1}{\color[rgb]{0,0,0}a_{i}\,\,\,\,}\left[{}_{t_{in}}d^{\alpha_{i}}_{t}y(t)+\Psi_{i}(f_{in},\alpha_{i},t_{abs},t_{in},t)\right]\\!\\!=b\,\,u(t),$
(28)

$f_{in}(t):\,\,\,[t_{abs}\,\,\,t_{in}]\rightarrow\mathbb{R}$ corresponds to
the history function of $y(t)$ based-initialization.

${\Psi=(\,\,\,\Psi_{1}(f_{in},\alpha_{1},t_{abs},t_{in},t)\,\,\,\,\,\,\Psi_{2}(f_{in},\alpha_{2},t_{abs},t_{in},t)\,\,\,\dots}$
${\,\,\,\Psi_{N}(f_{in},\alpha_{N},t_{abs},t_{in},t)\,\,\,)^{tr}}$ denotes the
vector of initialization functions and and are defined as in (9).

We denote the vectors $p$ and $\alpha$ as: $p=(a_{1}\ a_{2}\ \dots
a_{N},b)^{tr}$, $\alpha=(\alpha_{1}\ \alpha_{2}\ \dots\alpha_{N})^{tr}$ and
their estimates as $\hat{p}=(\hat{a}_{1}\ \hat{a}_{2}\
\dots\hat{a}_{N},\hat{b})^{tr}$ and $\hat{\alpha}=(\hat{\alpha}_{1}\
\hat{\alpha}_{2}\ \dots\hat{\alpha}_{N})^{tr}$ respectively. $(\cdot)^{tr}$
denotes the transpose of the row vector.

In this work, we are interested in solving the following estimation problem
(EP) that accounts for the initialized fractional-order derivative:

(EP) $\left\\{\begin{tabular}[]{l}{Given the output signal $y(t)$ and $u(t)$,
$t\in[t_{in},T]$,} {design a history function based-initialization,}\\\
{$f_{in}(t)$, $t\\!\\!\in[t_{abs},t_{in}]$ jointly with the simultaneous}
{estimation of the unknown parameters and}\\\ { the fractional orders, denoted
$(\hat{p},\hat{\alpha})$}\end{tabular}\right.$

##  4 Reformulation of the Estimation Problem in Discrete Space

In this section, we reformulate the initialized estimation problem in the
discrete space using the Grunwald-Letnikov derivative definition given in
Definition 4 and the history function based $GL$ given in Definition 5. In the
following, we present the different steps of the problem formulations:
For a given small enough sampling time $h$ such that
${h={(t-t_{abs})}\setminus{N_{1}}}$,
${N_{3}=\left({t_{in}-t_{abs}}\setminus{t-t_{abs}}\right)N_{1}}$ and
$K=N_{1}-N_{3}$. In the discrete space the system (28) can be written as
follows $\forall\,\,\,k=1,2,...K$

$y(t_{k})+{\sum}^{N}_{i=1}a_{i}\,\,\,\,\left[{}_{t_{in}}d^{\alpha_{i}}_{t}y(t_{k})+\Psi_{i}(f_{in},\alpha_{i},t_{abs},t_{in},t_{k})\right]\\!\\!=b\,\,u(t_{k}),$
(29)

where

$\begin{matrix}\Psi_{i}(f_{in},\alpha_{i},t_{in},t_{abs},t_{k})=\dfrac{1}{h^{\alpha_{i}}}{\sum}_{n=0}^{N_{3}-1}C_{n}^{\alpha_{i}}f_{in}[t_{k}-(N_{1}-1-n)h].\end{matrix}$
(30)

where $C_{j}^{\alpha_{i}}(j=0,1,\dots N_{1})$ are the binomial coefficients
recursively computed using the following formula:

$C_{0}^{\alpha_{i}}=1,\,\,\,\,\,\,\,C_{j}^{\alpha_{i}}=\left(1-\dfrac{1+\alpha}{j}\right)C_{j-1}^{\alpha_{i}}.$
(31)

Substituting (30) in (29) we obtain:

$y(t_{k})+{\sum}^{N}_{i=1}a_{i}\,\,\,\,\left[{}_{t_{in}}d^{\alpha_{i}}_{t}y(t_{k})+\dfrac{1}{h^{\alpha_{i}}}{\sum}_{n=0}^{N_{3}-1}C_{n}^{\alpha_{i}}f_{in}[t_{k}-(N_{1}-1-n)h]\right]\\!\\!=b\,\,u(t_{k}).$
(32)

Based on Definition 4 and substituting $d^{\alpha_{i}}_{t}y(t_{k})$ by its
formula described in (6) we obtain the following:

$y(t_{k})+{\sum}^{N}_{i=1}a_{i}\,\,\,\,\left[\dfrac{1}{h^{\alpha_{i}}}\,\displaystyle{\sum_{j=0}^{N_{1}}}C_{j}^{\alpha_{i}}y(t_{k-j})+\dfrac{1}{h^{\alpha_{i}}}{\sum}_{n=0}^{N_{3}-1}C_{n}^{\alpha_{i}}f_{in}[t_{k}-(N_{1}-1-n)h]\right]\\!\\!=b\,\,u(t_{k}).$
(33)

By rearranging (33) we obtain:

${\sum}^{N}_{i=1}a_{i}\,\,\,A_{i}(t_{k})+b\,\,u(t_{k})=y(t_{k})$ (34)

where

$A_{i}(t_{k})=-\dfrac{1}{h^{\alpha_{i}}}\left[\displaystyle{\sum_{j=0}^{N_{1}}}C_{j}^{\alpha_{i}}y(t_{k-j})+{\sum}_{n=0}^{N_{3}-1}C_{n}^{\alpha_{i}}f_{in}[t_{k}-(N_{1}-1-n)h]\right]$
(35)

Accordingly, the reformulated discrete initialized discretized estimation
problem (IDEP) is given in the following proposition.

######  Proposition 1.

For a given small enough sampling time $h$, the estimation problem (EP) can be
reformulated in the discrete space as follows:

$F(\alpha;N_{3})\cdot\,\,\,p=R,$ (36)

where $p$ is the vector of unknown parameters and $R\in\mathbb{R}^{K\times 1}$
such that $K=N_{1}-N_{3}$, the number of samples within $[t_{in},T]$,
$R_{k}=y(t_{k}),\,\,\,k=1,2,\dots,K$.
The function $F(\alpha;N_{3})\in\mathbb{R}^{K\times(N+1)}$ that is a non-
linear function with respect to the fractional differentiation orders and the
length of the initialization memory is given as follows: for
$\forall\,\,\,\,k=1:K$

$\displaystyle
F_{ki}(\alpha_{i};N_{3})=\left\\{\begin{array}[]{lll}A_{i}(t_{k})=-\dfrac{1}{h^{\alpha_{i}}}\left[\displaystyle{\sum_{j=0}^{N_{1}}}C_{j}^{\alpha_{i}}y(t_{k-j})+{\sum}_{n=0}^{N_{3}-1}C_{n}^{\alpha_{i}}f_{in}[t_{k}-(N_{1}-1-n)h]\right],\\\
\forall i=1:N\\\ u(t_{k}),\hskip 5.69046ptfor\,\,i=N+1,\end{array}\right.$
(40)

######  Remark 2.

It is worth noting that the reformulated problem (36) is expressed as a linear
combination of nonlinear functions. The unknown variables in the objective
function based on in (36), can be separated into two disjoint sets: linear
variables, $p$, and nonlinear variables, $(\alpha;N_{3})$.This class of
problems is known as the separable nonlinear least squares problem [20].

##  5 Estimation algorithm

Figure 4: A flowchart showing the two main steps of initialized estimation
algorithm for parameters and fractional differentiation orders of fractional-
order systems.

This section introduces the main two steps of the proposed method to solve
(EP). The first step designs a time-varying history function to initialize the
output signal of the FOS (27). The second step is an iterative hybrid approach
to solve the estimation problem of the unknown parameters and fractional
differentiation orders, (36). This step contains two stages. The first stage
solves the parameters’ estimation problem at each iteration using a linear
least square method by solving a linear system of equations; the second stage
solves a nonlinear system of equations using Newton’s method to find estimates
for the fractional differentiation orders. The main steps of the proposed
method which are explained in the following are depicted in Figure 4.

###  5.1 Step 1: Design of the time-varying history function

The key hypothesis in step 1 is to propose an initialization function,
$f_{in}(t)$, in the estimation process of the unknown parameters and
fractional differentiation orders. The choice of $f_{in}(t)$ depends on the
available measured output $y(t)$ where duplicates of $y(t)$ is considered.
This choice reduces the infinite-dimensional space of possible initialization
functions into a finite parametric space where the remaining degree of freedom
is the number of output cycles to be duplicated. A schematic picture of the
designed initialization framework (function) is depicted in Figure 4, Step 1.

###  5.2 Step 2: Two-stage Estimation Algorithm for Parameters and Fractional
Differentiation Order

Let the fractional derivatives of the output signal, $y(t)$, be initialized
with a time-varying function $f_{in}(t)$ as designed in the previous step. In
this study, we suppose that $N_{3}$ is known. Accordingly, in this situation
the set of parameters $p$ enters linearly while the set of parameters $\alpha$
non-linearly in the following predictor:

$\hat{y}(t\mid p,\alpha)=p^{tr}F^{tr}(t,\alpha;N_{3}).$ (41)

The identification criterion then becomes

$V_{K}(p,\alpha)={\sum}^{K}_{k=1}\mid
y(t_{k})-p^{tr}F^{tr}(t_{k},\alpha;N_{3})\mid^{2}=\mid
R-F(\alpha;N_{3})p\mid,$ (42)

where we introduced matrix notation, analogously to (36). Then, for a given
$\alpha$, this criterion is the least square criterion and minimized w.r.t.
$p$ by

$\hat{p}=\left[F^{tr}(\alpha;N_{3})F(\alpha;N_{3})\right]^{-1}F^{tr}(\alpha;N_{3})\,\,\,R$
(43)

By inserting (43) into the identification criterion (42) the non-linear
problem with repst to $\alpha$ can be defined as:

$\underset{\alpha}{\min}\mid
R-P(\alpha)R\mid^{2}=\underset{\alpha}{\min}\left|\left[I_{K\times
K}-P(\alpha)\right]R\right|^{2}\coloneqq\underset{\alpha}{\min}\mid
J(\alpha;N_{3})\mid^{2}$ (44)

where

$J(\alpha;N_{3})=\left[I_{K\times K}-P(\alpha)\right]R$ (45)

with $I_{K\times K}$ being the identity matrix of dimension $K$ and
$P(\alpha)$ is defined as:

$P(\alpha)=F(\alpha;N_{3})\left[F^{tr}(\alpha;N_{3})F(\alpha;N_{3})\right]^{-1}F^{tr}(\alpha;N_{3})$
(46)

######  Remark 3.

The estimate $\hat{p}$ is obtained by inserting the minimizing $\alpha$ into
(45). The matrix $P$ is a projection matrix: $P^{2}=P$. The described method
is called separable least squares since the least-square part has been
separated out, and the problem is reduced to a minimization problem of lower
dimension.

An estimate $\hat{\alpha}^{l}$ of the vector of fractional differentiation
orders are computed iteratively by solving the nonlinear problem:

$J(\alpha;N_{3})=0$ (47)

with $J\in\mathbb{R}^{K\times 1}$, which is defined in Eqt. (45), using the
following Newton update law:

$\hat{\alpha}^{l+1}=\hat{\alpha}^{l}-[J^{\prime}(\hat{\alpha}^{l};N_{3})]^{-1}J(\hat{\alpha}^{l};N_{3}),$
(48)

where $J^{\prime}(\alpha;N_{3})\in\mathbb{R}^{K\times N}$ is the Jacobian
matrix of the nonlinear function $J$.

###### Proof.

The system of equations (36) is solved in two stages. While in the first step
we use the linear least square method to estimate the parameters, in the
second step we use Newton iterative method to estimate the fractional
differentiation orders.
First stage: From (36) we have

$F({\alpha^{l}};N_{3})\cdot\,\,\,p=R,$ (49)

For any given $\alpha$, we denote $\hat{p}$ the solution of (49) for the
unknown vector of parameters $p$. Thus, for $\alpha=\hat{\alpha}^{l}$ the
system of equations(49) can be written as in (43).
Second stage: From (36) and the first stage we have

$F({\alpha};N_{3})\cdot\,\,\,\hat{p}=R,$ (50)

where $F({\alpha};N_{3})$ and $R$ are introduced in (40). We propose to use
the iterative first-order Newton method to solve(50), with a solution at
iteration $l$ denoted $\hat{\alpha}^{l}$. The Newton update law is given in
(48). At each iteration, (47) follow directly from (50).

The Jacobian matrix ($J^{\prime}(\alpha)$) needs to be evaluated at each
iteration $l$. Thus, we first substitute $\alpha$ by $\alpha^{l}$, and then we
use it in the updated law (48). ∎

The main steps of the proposed method are given in Algorithm 1.

Algorithm 1 Proposed two-stage estimation method

1:initialize $l=0$, give an initial guess to $\alpha^{l}$,

2:compute $\hat{p}$ using equation (43) ,

3:compute $J(\alpha^{l})$ using equations (45),

4:if $\parallel\dfrac{J(\alpha^{l})}{R}\parallel\prec\epsilon$
then222$\epsilon$ is a user-defined threshold that takes small values.

5:stop

6:else

7:update $\alpha^{l}$ using equations,

8:update the number of iterations $l=l+1$,

9:return to step 2

10:end if

##  6 Convergence of the algorithm

###  6.1 Characterization of the Jacobian matrix

In this part, we show the analytical derivation of the Jacobian matrix
$J^{\prime}(\alpha)$. The following lemmas are needed in the derivation of the
Jacobian matrix and the convergence of the proposed algorithm.

######  Lemma 1.

[21], Let $M\in R^{K\times N}$. $M^{tr}M$ is called the Gram matrix of $M$.
The Moore-Penrose inverse of $M$ is $M^{+}$ and we have:

$M^{+}=(M^{tr}M)^{+}M^{tr}=M^{tr}(M\,\,M^{tr})^{+}.$ (51)

######  Lemma 2.

Let $M(\alpha)$ be a defined on an interval $I$ whose values are nonzero
$n\times n$ matrices, differentiable at $\alpha^{*}$. Then the function
$M^{+}$ is differentiable at $\alpha^{*}$ if an only if $rank(M(\alpha))$ is
constant in some interval $|\alpha-\alpha^{*}|<\delta$. The derivative of the
$\dfrac{\partial M^{+}}{\partial\alpha}(\alpha^{*})$ is given by:

$\begin{array}[]{ll}\dfrac{\partial
M^{+}}{\partial\alpha}=-M^{+}\left(\dfrac{\partial
M}{\partial\alpha}\right)M^{+}+M^{+}(M^{+})^{tr}\left(\dfrac{\partial
M^{tr}}{\partial\alpha}\right)e+e\left(\dfrac{\partial
M^{tr}}{\partial\alpha}\right)(M^{+})^{tr}M^{+},\end{array}$ (52)

where $e=\left(I-MM^{+}\right)$ and $M,M^{tr},M^{+},\dfrac{\partial
M}{\partial\alpha}$ stand for $M(\alpha^{*})$, $M^{tr}(\alpha^{*})$,
$M^{+}(\alpha^{*})$, $\dfrac{\partial M}{\partial\alpha}(\alpha^{*})$
respectively.

Based on the formulation of $J(\alpha;N_{3})$ in (45) and (46) the elements of
the Jacobian matrix $J^{\prime}(\alpha;N_{3})$ can be computed as:
$\forall\,\,\,$ $i=1:K$ and, $j=1:N$:

$J^{\prime}_{ij}=\dfrac{\partial
J_{i}}{\partial\alpha_{j}}={\sum}_{l=1}^{K}\dfrac{\partial}{\partial\alpha_{j}}P_{il}(\alpha)R_{l}.$
(53)

Such that:

$P(\alpha)=F(\alpha)\,\,\,\overbrace{\left[\underbrace{F^{tr}(\alpha)F(\alpha)}_{Gram\,\,\,\,matrix}\right]^{-1}F^{tr}(\alpha)}^{F^{+}}.$
(54)

The matrix $F(\alpha)$ can be expressed as:

$F(\alpha)=\begin{pmatrix}&F_{11}(\alpha)&F_{12}(\alpha)&F_{13}(\alpha)&\dots&F_{1(N+1)}(\alpha)&\\\
&F_{21}(\alpha)&F_{22}(\alpha)&F_{23}(\alpha)&\dots&F_{2(N+1)}(\alpha)&\\\
&F_{31}(\alpha)&F_{32}(\alpha)&F_{33}(\alpha)&\dots&F_{3(N+1)}(\alpha)&\\\
&.&.&.&\dots&.&\\\ &.&.&.&\dots&.&\\\ &.&.&.&\dots&.&\\\
&F_{K1}(\alpha)&F_{K2}(\alpha)&F_{K3}(\alpha)&\dots&F_{K(N+1)}(\alpha)&\end{pmatrix}$
(55)

Based on Lemma. 1, we note:

$B=F^{+}=\left[F^{tr}(\alpha)F(\alpha)\right]^{-1}F^{tr}(\alpha)=\left[F^{tr}(\alpha)F(\alpha)\right]^{+}F^{tr}(\alpha),$
(56)

is the Moore-Penrose inverse matrix. Accordingly, (54) can be written as:

$P(\alpha)=F(\alpha)\,\,\,F^{+}(\alpha)=F(\alpha)\,\,\,B(\alpha),$ (57)

where $B(\alpha)$ is defined as:

$B(\alpha)=F^{+}(\alpha)=\begin{pmatrix}&B_{11}(\alpha)&B_{12}(\alpha)&B_{13}(\alpha)&\dots&B_{1K}(\alpha)&\\\
&B_{21}(\alpha)&B_{22}(\alpha)&B_{23}(\alpha)&\dots&B_{2K}(\alpha)&\\\
&B_{31}(\alpha)&B_{32}(\alpha)&B_{33}(\alpha)&\dots&B_{3K}(\alpha)&\\\
&.&.&.&\dots&.&\\\ &.&.&.&\dots&.&\\\ &.&.&.&\dots&.&\\\
&B_{(N+1)1}(\alpha)&B_{(N+1)2}(\alpha)&B_{(N+1)3}(\alpha)&\dots&B_{(N+1)K}(\alpha)&\\\
\end{pmatrix}$ (58)

The elements of the matrix $P(\alpha)$ are defined as follows:
$\forall\,\,\,i=1:K$ and $l=1:K$

$P_{il}(\alpha)={\sum}_{m=1}^{N+1}F_{im}(\alpha)\,\,\,B_{ml}(\alpha).$ (59)

and its derivative with respect to the fractional differentiation order
$\alpha$ can be evaluated as:

$\dfrac{\partial
P_{il}(\alpha)}{\partial\alpha_{j}}={\sum}_{m=1}^{N+1}F_{im}(\alpha)\dfrac{\partial
B_{ml}(\alpha)}{\partial\alpha_{j}}+{\sum}_{m=1}^{N+1}\dfrac{\partial
F_{im}(\alpha)}{\partial\alpha_{j}}B_{ml}(\alpha)$ (60)

where $\dfrac{\partial B_{ml}(\alpha)}{\partial\alpha_{j}}$ are the element of
the derivative of Moore-Penrose inverse matrix which can be computed using
Lemma. 2.

$\dfrac{\partial F_{im}(\alpha)}{\partial\alpha_{j}}$ is derived based on the
formulation of (40) as:

$\dfrac{\partial
F_{im}}{\partial\alpha_{j}}=\left\\{\begin{array}[]{ll}-\dfrac{\ln(h)}{h^{\alpha_{m}}}\left[{\sum}_{l=1}^{N_{1}}C_{l}^{\alpha_{m}}y(t_{i-l})+{\sum}_{n=1}^{N_{3}-1}C_{n}^{\alpha_{m}}f_{in}\left[t_{i}-(N_{1}-1-n)h\right]\right]-\\\
\\\ \dfrac{1}{h^{\alpha_{m}}}\left[{\sum}_{l=0}^{N_{1}}\dfrac{\partial
C_{l}^{\alpha_{m}}}{\partial\alpha_{m}}y(t_{i-l})+{\sum}_{n=0}^{N_{3}-1}\dfrac{\partial
C_{n}^{\alpha_{m}}}{\partial\alpha_{m}}f_{in}\left[t_{i}-(N_{1}-1-n)h\right]\right]&\hskip
28.45274ptj=m\\\ \\\ 0&\hskip 28.45274ptj\neq m\end{array}\right.$ (61)

where $\ln(.)$ is logarithmic function and the partial derivative of the
binomial coefficients with respect to $\alpha_{j}$ is computed iteratively as
follows:

$\displaystyle\dfrac{\partial
C_{n}^{\alpha_{j}}}{\partial\alpha_{j}}=\left\\{\begin{array}[]{lll}0\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,n=0\\\
\\\
\dfrac{-1}{n}C_{n-1}^{\alpha_{j}}+\left(1-\dfrac{1+\alpha_{j}}{n}\right)\dfrac{\partial
C_{n-1}^{\alpha_{i}}}{\partial\alpha_{j}}\,\,\,\,\,\,\,\,\,n\geq
1\end{array}\right.$ (65)

Substituting (60) in (53) the element of the Jacobian matrix,
$J^{\prime}(\alpha)$ can be computed using the following expression:

$J^{\prime}_{ij}(\alpha)={\sum}_{l=1}^{K}\,\,\,{\sum}_{m=1}^{N+1}F_{im}(\alpha)\dfrac{\partial
B_{ml}}{\partial\alpha_{j}}(\alpha)+{\sum}_{l=1}^{K}\dfrac{\partial
F_{ij}}{\partial\alpha_{j}}B_{jl}$ (66)

The matrix $J^{\prime}$ can be expressed as:

$J^{\prime}(\alpha)=\begin{pmatrix}&\dfrac{\partial
J_{1}}{\partial\alpha_{1}}&\dfrac{\partial
J_{1}}{\partial\alpha_{2}}&\dfrac{\partial
J_{1}}{\partial\alpha_{3}}&\dots&\dfrac{\partial
J_{1}}{\partial\alpha_{N+1}}&\\\ &.&.&.&\dots&.&\\\ &.&.&.&\dots&.&\\\
&.&.&.&\dots&.&\\\ &\dfrac{\partial J_{K}}{\partial\alpha_{1}}&\dfrac{\partial
J_{K}}{\partial\alpha_{2}}&\dfrac{\partial
J_{K}}{\partial\alpha_{3}}&\dots&\dfrac{\partial
J_{K}}{\partial\alpha_{N+1}}&\end{pmatrix}$ (67)

###  6.2 Proof of convergence

The iterative character of the above algorithm results from the Newton method
which is used to estimate the fractional differentiation order, $\alpha$. The
vector of parameters $p$ is estimated by solving the least square problem.
However, we note that this estimate depends on the vector of fractional orders
$\alpha$.
Thus, we use a standard convergence result of Newton’s method to provide
sufficient conditions to ensure the local convergence of the proposed two-
stage algorithm. Such a result can be found in Chapter 5 of [22]. We also make
use of two results on Lipschitz functions stating that the product and
quotient (denominator not zero) of two bounded Lipschitz functions are
Lipschitz. The following assumptions are needed to prove the convergence
result stated in Theorem 1.

######  Assumption 1.

The input $u(t)$, output $y(t)$ and the initialization function $f_{in}(t)$
are bounded functions in $[0,T]$.

######  Assumption 2.

A solution $\alpha^{*}$ to the system $J(\alpha)=0$ exists.

######  Assumption 3.

The Jacobian matrix $J^{\prime}(\alpha^{*})$ is nonsingular.

######  Theorem 1.

Assume that Assumptions $\ref{assum1}-\ref{assum3}$ hold. Then, the proposed
proposed two-stage algorithm converges q-quadratically.

Proof. Based on Theorem 5.1.2, page 71 in [22], and given Assumptions 2 and 3,
it remains to prove that the Jacobian matrix is Lipschitz continuous for
$\alpha\in\wp=(n_{i-1},n_{i})$
333$\wp=\prod^{N}_{i=1}\wp_{i}=\\{(a_{1},a_{2}\dots,a_{N})|\forall
i=1,2,\dots,N,a_{i}\in\wp_{i}\\}$ is the generalized Cartesian product of $N$
sets $\wp_{1},\wp_{2}\dots\wp_{N}$. in order to conclude about the convergence
of the proposed two-stage algorithm.

The entries of the Jacobian matrix are given by (66) where the terms $F_{im}$,
$\dfrac{\partial F_{im}(\alpha)}{\partial\alpha_{j}}$, $B_{jl}$ and
$\dfrac{\partial B_{ml}(\alpha_{j})}{\partial\alpha_{j}}$ are characterized in
(40), (61), (56) and Lemma 2, respectively.

The norm 1 of the term $J^{\prime}(\alpha^{1})-J^{\prime}(\alpha^{2})$ is
given as follows:

$\displaystyle\begin{array}[]{l}\parallel
J^{\prime}(\alpha^{2})-J^{\prime}(\alpha^{2})\parallel_{1}=\underset{1\leq
j\leq N}{\max}{\sum}_{i=1}^{K}\mid
J^{\prime}_{ij}(\alpha^{2})-J^{\prime}_{ij}(\alpha^{2})\mid\\\ \\\ \hskip
99.58464pt=\underset{j}{\max}{\sum}_{i=1}^{K}\Bigg{|}\left({\sum}_{l=1}^{K}\,\,\,{\sum}_{m=1}^{N+1}F_{im}(\alpha^{1})\dfrac{\partial
B_{ml}}{\partial\alpha_{j}}(\alpha^{1})+{\sum}_{l=1}^{K}\dfrac{\partial
F_{ij}}{\partial\alpha_{j}}B_{jl}\right)\\\ \\\ \hskip
113.81102pt-{\sum}_{l=1}^{K}\,\,\,{\sum}_{m=1}^{N+1}F_{im}(\alpha^{2})\dfrac{\partial
B_{ml}}{\partial\alpha_{j}}(\alpha^{2})+{\sum}_{l=1}^{K}\dfrac{\partial
F_{ij}}{\partial\alpha_{j}}(\alpha^{2})B_{jl}\Bigg{|}\\\ \\\ \hskip
99.58464pt={\sum}_{i=1}^{K}\Bigg{|}{\sum}_{l=1}^{K}\,\,\,{\sum}_{m=1}^{N+1}\left\\{F_{im}(\alpha^{1})\dfrac{\partial
B_{ml}}{\partial\alpha_{j^{*}}}(\alpha^{1})-F_{im}(\alpha^{2})\dfrac{\partial
B_{ml}}{\partial\alpha_{j^{*}}}(\alpha^{2})\right\\}+\\\ \\\ \hskip
113.81102pt\left({\sum}_{l=1}^{K}\left\\{\dfrac{\partial
F_{ij^{*}}}{\partial\alpha_{j^{*}}}(\alpha^{1})B_{j^{*}l}(\alpha^{1})-\dfrac{\partial
F_{ij^{*}}}{\partial\alpha_{j^{*}}}(\alpha^{2})B_{j^{*}l}(\alpha^{2})\right\\}\right)\Bigg{|}\end{array}$
(77)

where $j^{*}$ is the column index corresponding to the maximum of the sum of
the row components. Then, using the triangle inequality, (77) can be written
as follows:

$\begin{array}[]{ll}\parallel
J^{\prime}(\alpha^{2})-J^{\prime}(\alpha^{2})\parallel_{1}\,\,\leqslant\\\ \\\
{\sum}_{i=1}^{K}\,\,\,{\sum}_{l=1}^{K}\,\,\,{\sum}_{m=1}^{N+1}\left|F_{im}(\alpha^{1})\dfrac{\partial
B_{ml}}{\partial\alpha_{j^{*}}}(\alpha^{1})-F_{im}(\alpha^{2})\dfrac{\partial
B_{ml}}{\partial\alpha_{j^{*}}}(\alpha^{2})\right|\\\ \\\ \hskip
8.5359pt+{\sum}_{i=1}^{K}\,\,\,{\sum}_{l=1}^{K}\left|\dfrac{\partial
F_{ij^{*}}}{\partial\alpha_{j^{*}}}(\alpha^{1})B_{j^{*}l}(\alpha^{1})-\dfrac{\partial
F_{ij^{*}}}{\partial\alpha_{j^{*}}}(\alpha^{2})B_{j^{*}l}(\alpha^{2})\right|\end{array}$
(78)

Given the right hand side of inequality (78) and using properties on Lipschitz
functions, if $F_{im}(\alpha)$ and $\dfrac{\partial
F_{ij}}{\partial\alpha_{j}}(\alpha)$ are Lipschitz continuous and
$B_{jl}(\alpha)$ and $\dfrac{\partial B_{ml}}{\partial\alpha_{j}}(\alpha)$ are
bounded then we can conclude about the Lipschitz continuity of
$J^{\prime}(\alpha)$.

Given that $C_{l}^{\alpha}$ and $\dfrac{\partial
C_{l}^{\alpha}}{\partial\alpha}$ is high-order polynomial of $\alpha$ then
they are Lipschitz continuous. In addition,
$\left(\dfrac{1}{h^{\alpha}}\right)$ is bounded on bounded local interval of
$\alpha$. Accordingly, and taking into account  Assumption 1 we can insure the
Lipschitz continuity of $F_{im}(\alpha)$ and $\dfrac{\partial
F_{ij}}{\partial\alpha_{j}}(\alpha)$.

With respect to the boundedness of $B_{jl}(\alpha)$ and $\dfrac{\partial
B_{ml}}{\partial\alpha_{j}}(\alpha)$, based on the proof in [] that confirm
that if $F(\alpha)$ has a constant rank for
$\alpha^{*}-\delta<\alpha\leq\alpha^{*}+\delta$ then its Moore-Penrose inverse
$B$ is differentiable at $\alpha^{*}$ and its derivative is computed using
Lemma 2. Accordingly the of $B_{jl}$ and its derivative follows from this
condition.

##  7 Numerical Results

The performance of the proposed method is illustrated through different
numerical examples. Furthermore, potential applications of the algorithm are
presented, which consists of estimating parameters and fractional
differentiation orders of 1) fractional-order two-element arterial Windkessel
(F-WK2) model and 2) fractional-order neurovascular model. In the following
analysis of the results, we consider the percent relative error $Re(\%)$ as a
measure of estimation accuracy. In the case of parameter $p$ estimation, the
$Re_{p}(\%)$ is defined as:

$Re_{p}=\dfrac{|\hat{p}-p|}{|p|}\times 100,$ (79)

where $\hat{p}$ denotes the estimate of $p$. In the case of signal estimation
$y$, the $Re_{y}(\%)$ is defined as:

$Re_{y}=\dfrac{||\hat{y}-y||_{2}}{||y||_{2}}\times 100,$ (80)

where $\hat{y}$ denotes the estimate of the real measurement $y$.

###  7.1 Example 1

Figure 5: Estimated output signal along with the real one using random input
signal. Here we use $N_{0}=15$ and $N_{C}=10$. The sampling rate $h=0.01$ has
been used in this simulation.

We consider the following fractional-order system:

$y(t)+a_{1}\,\,\,{}_{t_{in}}D_{t}^{\alpha}y(t)=b\,\,\,u(t)\,\,\,\,\,\,t\in[0,\,\,\,0.84]$
(81)

Here we consider two numerical cases based on the nature of the input signal
$u(t)$, namely: Random input signal and periodic pulse train input signal.

####  7.1.1 Random input signal

In this example, we consider the system described in (81) with a random input
signal whose amplitude varies between $0$ and $100$. The parameters and
fractional differentiation order are $a_{1}=1$, $b=0.5$ and $\alpha=0.7$. A
sensitivity analysis of the proposed algorithm with respect to the length of
the time-varying history function-based output fractional-order derivative
initialization, $N_{3}$, and the length of the original signal was performed.
In addition, a sensitivity analysis with respect to the initial guess of
Newton’s method (the initial guess of the fractional differentiation order
$\alpha^{0}$) is performed.

Figure 6: Contour plot of the relative error (%) with respect to the number of
cycles based time-varying history function, $N_{C}$ and the number of original
cycles based for estimation, $N_{0}$ for: parameter ($a_{1}$), parameter
($b$), the fractional differentiation order ($\alpha$) and the reconstructed
output signal ($y$). The sampling rate $h=0.01$ has been used in this
simulation. Figure 7: Relative error (%) with respect to the initial guess of
the fractional differentiation order $\alpha_{0}$ for: parameter ($a_{1}$),
parameter ($b$), the fractional differentiation order ($\alpha$) and the
reconstructed output signal ($y$).The sampling rate $h=0.01$ has been used in
this simulation.

We note $N_{C}=N_{3}/L$ and $N_{0}=K/L$ where $L$ corresponds to the length of
the original output signal and $K$ is the total number of samples used for
estimation, following the previous notations, it is defined as
$K=N_{1}-N_{3}$.
Figure 5 shows the input signal and the estimated output signal along with the
original one for $N_{0}=15$ and $N_{C}=10$. In this case, the relative errors
of the parameters estimate $a_{1}$ and $b$ were $11.77\%$ and $0.13\%$
respectively, whereas, for the fractional differentiation order estimate, it
was $3.93\%$ and the reconstructed output around $1.22\%$. It is clear from
these results that we have a good estimation performance for the output
signal. Although the estimation error of the parameters and the fractional
differentiation order is not the best, we still have a good reconstruction of
the output. The reason is that in this case, we cannot insure the global
identifiability of the system (81).
Figure 6 shows a contour plot of the relative errors (%) with respect to the
number of cycles based time-varying history function, $N_{C}$, and the number
of original cycles based for estimation, $N_{0}$ for the parameters,
fractional differentiation order and the reconstructed output. Overall, from
these subplots it is clear that increasing both $N_{C}$ and $N_{0}$ leads to a
decrease in the relative error and so better estimation for the parameters,
fractional differentiation order, and the reconstructed output. For a fixed
$N_{C}$ we noticed that as $N_{0}$ increases the relative error decreases.
With regards to the unknown parameters and $\alpha$, there are optimal values
of ($N_{C}^{*},N_{0}^{*}$) that lead to the minimum error. Any combination
larger than ($N_{C}^{*},N_{0}^{*}$) might lead to a small increase in the
error. With respect to the output relative error, we noticed that increasing
more $N_{0}$ leads to better performances. For small values of $N_{0}=1,2\dots
5$ it is clear that there is the optimal combination ($N_{C}^{*},N_{0}^{*}$)
that leads to the best performance.
Figure 7 summarizes the results of the relative errors of the estimated
unknowns when the initial guess of Newton’s method is varying. In this
simulation we varied $\alpha^{0}$ from $0.05$ to $2$ with a step of $0.05$. It
is clear that the proposed algorithm is not very sensitive to initial guesses
in this case.

Figure 8: Estimated output signal along with the real one using random input
signal. Here we use $N_{0}=15$ and $N_{C}=10$. The sampling rate $h=0.01$ has
been used in this simulation. Figure 9: Contour plot of the relative error (%)
with respect to the number of cycles based time-varying history function,
$N_{C}$ and the number of original cycles based for estimation, $N_{0}$ for:
parameter ($a_{1}$), parameter ($b$), the fractional differentiation order
($\alpha$) and the reconstructed output signal ($y$). The sampling rate
$h=0.01$ has been used in this simulation. Figure 10: Relative error (%) with
respect to the initial guess of the fractional differentiation order
$\alpha_{0}$ for: parameter ($a_{1}$), parameter ($b$), the fractional
differentiation order ($\alpha$) and the reconstructed output signal ($y$).
The sampling rate $h=0.01$ has been used in this simulation.

####  7.1.2 Periodic pulse train input signal

In this example we consider the system described in (81) where the parameters
and fractional differentiation order are $a_{1}=1$, $b=0.5$ and $\alpha=0.7$.
The input $u(t)$ is taken as a periodic pulse train of duration
$T=0.84\,\,\,[s]$, an amplitude $A=350$ and a duty cycle ($d_{C}=0.27$). The
time step used in this example is $h=0.01$. Similar to the previous example, a
sensitivity analysis of the proposed algorithm with respect to the length of
the time-varying history function-based output fractional-order derivative
initialization, the length of the original signal and the initial guess of
Newton’s method is performed. Figure 8 summarizes the reconstruction results
of the output of the fractional-order system (81) using a periodic pulse train
signal as an output. In this simulation we used ($N_{0}=15;N_{C}=10$). The
relative errors of this reconstruction is $1.19\%$, for the unknown parameter
$a_{1}$ is $2.42\%$ and $b$ is $1.56\%$ and the fractional-differentiation
order is $0.35\%$.
In figure 9, we summarize the results of the relative errors of the estimated
unknowns and the reconstructed output for a different combination of
($N_{0};N_{C}$).
Figure 10 summarizes the results of the relative errors of the estimated
unknowns when the initial guess of the fractional differentiation order,
$\alpha^{0}$, is varying from $0.4$ to $1.25$ with a step $0.05$.

####  7.1.3 Joint estimation of the parameters and the fractional
differentiation order of F-WK2

A potential application of the proposed algorithm is to estimate the
hemodynamic parameters and the differentiation order of the fractional-order
two-element arterial Windkessel (F-WK2). As shown in Fig. 11, F-WK2 represents
the heart as a current source that pumps blood into the arterial system, which
is lumped into a single fractional-order capacitor and a single resistance
[23, 24]. The fractional-order capacitor, characterized by two parameters
($C_{\alpha},\alpha$), represents the apparent arterial compliance by taking
into account the viscoelastic properties of the vascular network via the
fractional-order factor [25, 26]. The resistance $R_{p}$ accounts for the
total peripheral arterial resistance. The F-WK2 system is given as follows:

$P_{a}(t)+\tau_{\alpha}\,\,\,{}_{t_{in}}D_{t}^{\alpha}P_{a}(t)=R_{p}\,\,\,Q_{a}(t),\,\,\,t\in[0,T_{C}]$
(82)

where $\small\tau_{\alpha}=C_{\alpha}\cdot R_{p}$, $\small 0\leq\alpha\leq 1$
and $\small T_{C}$ is the cardiac period. $\small P_{a}$ corresponds to the
aortic blood pressure, and $\small Q_{a}$ denotes the aortic blood flow.

Figure 11: Fractional-order tow-element arterial Windkessel circuit model.
$C_{\alpha}$ represents the arterial compliance and $R_{p}$ represents the
peripheral resistance. Figure 12: Estimated aortic blood pressure using F-WK2
along with the real one. Here we use $N_{0}=10$ and $N_{C}=25$. The relative
errors of the reconstructed output is $5.22\%$, $\tau_{\alpha}=1.15$ and
$R_{p}=1.13$. The sampling rate $h=0.01$ has been used in this simulation.

To fully identify F-WK2, the hemodynamic parameters and the fractional
differentiation order have to be estimated from measured flow and pressure
waves from accessible arterial locations. This is also known as the
hemodynamic inverse problem [27]. The proposed algorithm has been applied to
an in-silico data generated from one-dimensional model [28]. Here, the aortic
blood pressure is considered the output of the fractional-order system, while
the aortic blood flow is considered the input. To initialize the fractional-
order derivative of the output ($D_{t}^{\alpha}P_{a}(t)$), we used a certain
number of cycles, $N_{C}$ of the same output signal as a time-varying history
function. We denote by $N_{0}$ the number of original cardiac cycles used for
estimation. A sensitivity analysis of the applied algorithm with respect to
$N_{C}$ and ($N_{0}$) is performed.
Fig. 12 summarizes the reconstruction results of $10$ cardiac cycles of the
aortic blood pressure using a real input and the estimated parameters and the
fractional differentiation order. This reconstruction used $25$ cycles to
initialize the fractional-order derivative of the output and the relative
error was around $5.22\%$. The estimate of the constant time is
$\hat{\tau}_{\alpha}=1.15$ and the peripheral resistance $\hat{R}_{p}=1.13$.
Figure 14 summarizes the sensitivity analysis of the relative error of the
reconstructed aortic blood pressure (output of the system (82)) with respect
$N_{C}$ and $N_{0}$. It is clear that as we increase $N_{C}$ the relative
error of the estimation decreases. Figure 14 summarizes the results of the
relative errors of the estimated output the initial guess is varying from
$0.35$ to $1.2$ with a step of $0.05$

Figure 13: Relative estimation error vs. number of cycles based time-varying
history function, $N_{C}$ for one cycle output aortic blood pressure signal
($N_{0}=1$). The sampling rate $h=0.01$ has been used in this simulation.

Figure 14: Relative error (%) with respect to the initial guess of the
fractional differentiation order $\alpha_{0}$ for the reconstructed output
signal, the aortic blood pressure. The sampling rate $h=0.01$ has been used in
this simulation.

###  7.2 Example 2

We consider the following example:

$y(t)+a_{1}\,D^{\alpha_{1}}_{t}y(t)+a_{2}\,D^{\alpha_{2}}_{t}y(t)=u(t),\,\,\,t\in[0,10]$
(83)

where the parameters and fractional differentiation orders are $a_{1}=3$,
$a_{2}=2$ and $\alpha_{1}=1.5$, $\alpha_{2}=0.5$. Here we consider two
numerical cases based on the nature of the input signal $u(t)$, namely: The
$sinc$ function as the input signal and the square input signal.

####  7.2.1 $Sinc$ function as an input signal

Figure 15: Estimated output signal along with the real one using random input
signal. Here we use $N_{0}=3$ and $N_{C}=10$. The sampling rate $h=0.1$ has
been used in this simulation. Figure 16: Contour plot of the percent relative
error (%) with respect to $N_{C}$ and $N_{0}$ for parameters, fractional
differentiation orders and the reconstructed output. The sampling rate $h=0.1$
has been used in this simulation.

In this case, we consider the following profile for the input $u(t)$ in the
system (83):

$u(t)=10\operatorname{sinc}(2\pi t)$ (84)

Figure 15 shows the reconstruction result of the estimated output using
$N_{C}=10$ and $N_{0}=3$. The relative errors of the unknown parameters
($a_{1},a_{2},\alpha_{1},\alpha_{2}$) are ($1.45\%,1.60\%,0.33\%,3.61\%$)
respectively and for the reconstructed output is $0.88\%$.
In figures 16, we summarize the results of the relative errors of the
estimated unknowns parameters, fractional differentiation orders, and the
reconstructed output for a different combination of ($N_{0};N_{C}$). From
these results, it is clear that as the $N_{C}$ increases the estimation
accuracy increases. It is worth noting that for a small value of $N_{C}$
increasing $N_{0}$ enhances the estimation accuracy.

####  7.2.2 Square input signal

Figure 17: Estimated output signal along with the real one using random input
signal. Here we use $N_{0}=3$ and $N_{C}=15$. Figure 18: Contour plot of the
relative error (%) with respect to $N_{C}$ and $N_{0}$ for parameters,
fractional differentiation orders and the reconstructed output signal. The
sampling rate $h=0.01$ has been used in this simulation.

In this example the input $u(t)$ of the fractional-order system (83) is taken
as a square function activated at $t=2$ and deactivated at $t=7$. Figure 17
shows the reconstruction result of the estimated output using $N_{C}=N_{0}=3$.

The relative errors of the unknown parameters
($a_{1},a_{2},\alpha_{1},\alpha_{2}$) are (4.19%, 4.17%, 1.07%, 5.84%)
respectively and for the reconstructed output is $0.86\%$. It is clear in this
example that although the percent relative error of the reconstructed output
is very low, the relative errors of the unknown parameters and fractional
differentiation estimates are higher. This can be explained by the fact that
fractional-order systems with more than one fractional differentiation order
are not structurally globally identifiable [29]. Figure 18 shows the percent
relative errors of the estimation of the parameters ($a_{1};a_{2}$), the
fractional differentiation orders ($\alpha_{1};\alpha_{2}$) and the
reconstructed output $y(t)$. Based on these results, overall, we have a good
reconstruction. In all cases, the percent relative error does not exceed
$3\%$. The best performance is achieved with large values of $N_{C}$ and
$N_{0}$. Concerning the parameters and fractional differentiation orders, we
obtained the same aspect. However, the error is larger than the ones obtained
in the reconstruction of the output.

####  7.2.3 Joint estimation of the parameters and the fractional
differentiation order for a fractional neurovascular mode

Figure 19: Estimated output signal along with the real one using random input
signal. Here we use $N_{0}=3$ and $N_{C}=10$.

A potential application of the proposed method is to estimate the neural
activity and the fractional differentiation orders from the Cerebral Blood
Flow (CBF) measurements. In [30], a fractional model for the neurovascular
coupling in the brain has shown a better characterization of the cerebral
hemodynamic response. This model relates the neural activity, considered as
input, to the CBF considered as output. To fully identify the fractional
neurovascular model, the input and the fractional differentiation orders have
to be estimated from CBF data. The nature of the neural activity can be of two
types: the block design paradigm, which corresponds to tasks lasting for
extended periods of time, and the event-related paradigm, which corresponds to
tasks lasting for shorter periods. This section focuses on estimating the
event-related type of neural activity usually modeled by Gaussian functions.
The neurovascular model is given as follows:

$D^{\alpha_{1}}_{t}f(t)+k\,D^{\alpha_{2}}_{t}f(t)+\gamma
f(t)=u(t),\,\,\,t\in[0,15]$ (85)

where $1\leq\alpha_{1}\leq 2$, $0\leq\alpha_{2}\leq 1$ are the fractional
differentiation orders, $f(t)$ is the temporal CBF and $u(t)$ is the temporal
neuronal activity. The constants $k$ and $\gamma$ are the flow signal decay
rate and the flow-dependant elimination constant, respectively. Their values
are given in as: $k=0.65$, $\gamma=0.41$.

Figure 20: Contour plot of the relative error (%) with respect to $N_{C}$ and
$N_{0}$ for parameters, fractional differentiation orders and the
reconstructed output signal. The sampling rate $h=0.1$ has been used in this
simulation.

The tested input profile is Gaussian which can be expressed as:

$u(t)=\exp\left[-(t-5)^{2}\right]$ (86)

We use a synthetic data set which have been generated from discretizing the
neurovascular model (85) and fixing ${\alpha_{1}=1.7,\alpha_{2}=0.6}$.

Figure 19 shows the reconstruction result of the estimated output using
$N_{C}=10$ and $N_{0}=3$. The relative errors of the unknown parameters
($k$,$\gamma$,$\alpha_{1},\alpha_{2}$) are ($1.32\%,1.64\%,0.62\%,1.67\%$)
respectively and for the reconstructed output is $0.57\%$.
In figures 20, we summarize the results of the relative errors of the
estimated unknown parameters, fractional differentiation orders, and the
reconstructed output for a different combination of ($N_{0};N_{C}$). From
these results, it is clear that as the $N_{C}$ increases the estimation
accuracy increases. It is worth noting that for a small value of $N_{C}$
increasing $N_{0}$ enhances the estimation accuracy.

##  8 Conclusion

The concept of initialized fractional calculus is of great pertinence for the
uniqueness of solutions to fractional forward problems. This report
investigates their significance for solving inverse problems. We proposed an
initialized estimation algorithm to jointly estimate the parameters and the
fractional differentiation orders while designing an output-dependent history
function based-initialization. After reformulating the estimation problem in
the discrete space, we proposed an algorithm of two steps that study the
effect of the length of the initialization function while estimating the
unknown parameters. We applied the pipeline to different numerical examples.
Furthermore, potential applications of the proposed algorithm are presented,
which consists of joint estimation of parameters and fractional
differentiation orders of a fractional-order arterial Windkessel and
neurovascular models. In future work, we will propose a multi-stage estimation
algorithm based on modulating functions that are known to be robust to noise.

## References

  * [1] R. Gorenflo and F. Mainardi, “Fractional calculus,” in _Fractals and fractional calculus in continuum mechanics_. Springer, 1997, pp. 223–276.
  * [2] R. L. Magin, “Fractional calculus models of complex dynamics in biological tissues,” _Computers & Mathematics with Applications_, vol. 59, no. 5, pp. 1586–1593, 2010.
  * [3] R. Matusu, “Application of fractional order calculus to control theory,” _International journal of mathematical models and methods in applied sciences_ , vol. 5, no. 7, pp. 1162–1169, 2011.
  * [4] V. E. Tarasov, “On history of mathematical economics: Application of fractional calculus,” _Mathematics_ , vol. 7, no. 6, p. 509, 2019.
  * [5] I. Podlubny, “Fractional-order systems and fractional-order controllers,” _Institute of Experimental Physics, Slovak Academy of Sciences, Kosice_ , vol. 12, no. 3, pp. 1–18, 1994.
  * [6] C. F. Lorenzo and T. T. Hartley, “Initialization of fractional-order operators and fractional differential equations,” _Journal of computational and nonlinear dynamics_ , vol. 3, no. 2, 2008.
  * [7] M. A. Bahloul, Z. Belkhatir, and T.-M. L. Kirati, “Initialization of fractional order systems for the joint estimation of parameters and fractional differentiation orders,” in _2022 American Control Conference (ACC)_. IEEE, 2022, pp. 3943–3948.
  * [8] G. B. H. Frej, R. Malti, M. Aoun, and T. Raissi, “Fractional interval observers and initialization of fractional systems,” _Communications in Nonlinear Science and Numerical Simulation_ , vol. 82, p. 105030, 2020.
  * [9] H. Ye and R. Huang, “Initial value problem for nonlinear fractional differential equations with sequential fractional derivative,” _Advances in Difference Equations_ , vol. 2015, no. 1, pp. 1–13, 2015.
  * [10] C. F. Lorenzo and T. T. Hartley, “Initialization in fractional order systems,” in _2001 European Control Conference (ECC)_. IEEE, 2001, pp. 1471–1476.
  * [11] ——, “Initialization, conceptualization, and application in the generalized (fractional) calculus,” _Critical Reviews™ in Biomedical Engineering_ , vol. 35, no. 6, 2007.
  * [12] B. Du, Y. Wei, S. Liang, and Y. Wang, “Estimation of exact initial states of fractional order systems,” _Nonlinear Dynamics_ , vol. 86, no. 3, pp. 2061–2070, 2016.
  * [13] D. Liu, T. M. Laleg-Kirati, O. Gibaru, and W. Perruquetti, “Identification of fractional order systems using modulating functions method,” in _2013 IEEE American Control Conference_. IEEE, 2013.
  * [14] Z. Belkhatir and T. M. Laleg-Kirati, “Joint estimation of the fractional differentiation orders and the unknown input for linear fractional non-commensurate system,” in _2015 IEEE Conference on Control Applications (CCA)_. IEEE, 2015, pp. 388–393.
  * [15] ——, “Parameters and fractional differentiation orders estimation for linear continuous-time non-commensurate fractional order systems,” _Systems & Control Letters_, vol. 115, pp. 26–33, 2018.
  * [16] M. A. Bahloul, M. Benencase, Z. Belkhatir, and T.-M. L. Kirati, “Finite-time simultaneous estimation of aortic blood flow and differentiation order for fractional-order arterial windkessel model calibration,” _IFAC-PapersOnLine_ , vol. 54, no. 15, pp. 538–543, 2021.
  * [17] S. Das, “Initialized differintegrals and generalized calculus,” in _Functional Fractional Calculus_. Springer, 2011, pp. 271–322.
  * [18] K. Diethelm, _The analysis of fractional differential equations: An application-oriented exposition using differential operators of Caputo type_. Springer Science & Business Media, 2010.
  * [19] M. Ö. Efe, “Fractional order systems in industrial automation—a survey,” _IEEE Transactions on Industrial Informatics_ , vol. 7, no. 4, pp. 582–591, 2011.
  * [20] C. F. Borges, “A full-newton approach to separable nonlinear least squares problems and its application to discrete least squares rational approximation,” _Electronic Transactions on Numerical Analysis_ , vol. 35, pp. 57–68, 2009.
  * [21] T. Kurmayya and K. Sivakumar, “Moore-penrose inverse of a gram matrix and its nonnegativity,” _Journal of optimization theory and applications_ , vol. 139, no. 1, pp. 201–207, 2008.
  * [22] C. T. Kelley, _Iterative methods for linear and nonlinear equations_. SIAM, 1995.
  * [23] M. A. Bahloul and T.-M. Laleg-Kirati, “Assessment of fractional-order arterial windkessel as a model of aortic input impedance,” _IEEE Open Journal of Engineering in Medicine and Biology_ , vol. 1, pp. 123–132, 2020.
  * [24] T. M. Laleg and M. A. Bahloul, “Mathematical biomarker for arterial viscoelasticity assessment,” Nov. 18 2021, uS Patent App. 17/284,072.
  * [25] M. A. Bahloul, Y. Aboelkassem, and T.-M. Laleg-Kirati, “Human hypertension blood flow model using fractional calculus,” _Frontiers in Physiology_ , vol. 13, p. 838593, 2022.
  * [26] ——, “Towards characterization of the complex and frequency-dependent arterial compliance based on fractional-order capacitor,” in _2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_. IEEE, 2021, pp. 5559–5565.
  * [27] S. Bernhard, K. Al Zoukra, and C. Schtte, “From non-invasive hemodynamic measurements towards patient-specific cardiovascular diagnosis,” in _Quality Assurance in Healthcare Service Delivery, Nursing and Personalized Medicine: Technologies and Processes_. IGI Global, 2012, pp. 1–25.
  * [28] P. H. Charlton, J. Mariscal Harana, S. Vennin, Y. Li, P. Chowienczyk, and J. Alastruey, “Modeling arterial pulse waves in healthy aging: a database for in silico evaluation of hemodynamics and pulse wave indexes,” _American Journal of Physiology-Heart and Circulatory Physiology_ , vol. 317, no. 5, pp. H1062–H1085, 2019.
  * [29] S. Alavi, A. Mahdi, P. E. Jacob, S. J. Payne, and D. A. Howey, “Structural identifiability analysis of fractional order models with applications in battery systems,” _arXiv preprint arXiv:1511.01402_ , 2015.
  * [30] Z. Belkhatir and T. M. Laleg-Kirati, “Fractional dynamical model for neurovascular coupling,” in _2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society_. IEEE, 2014, pp. 4916–4919.