# Output error MISO system identification using fractional models

Abir Mayoufi, Stephane Victor, Manel Chetoui, Rachid Malti, Mohamed Aoun

► To cite this version:

Abir Mayoufi, Stephane Victor, Manel Chetoui, Rachid Malti, Mohamed Aoun. Output error MISO system identification using fractional models. Fractional Calculus and Applied Analysis, 2021, 24 (5), pp.1601-1618. 10.1515/fca-2021-0067. hal-03548488

HAL Id: hal-03548488

https://hal.science/hal-03548488

Submitted on 1 Feb 2024

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.
# Output error MISO system identification using fractional models

Abir Mayoufi, Stéphane Victor, Manel Chetoui, Rachid Malti, Mohamed Aoun

► To cite this version:

Abir Mayoufi, Stéphane Victor, Manel Chetoui, Rachid Malti, Mohamed Aoun. Output error MISO system identification using fractional models. Fractional Calculus and Applied Analysis, 2021, 24 (5), pp.1601-1618. 10.1515/fca-2021-0067. hal-03548488

HAL Id: hal-03548488

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

Submitted on 1 Feb 2024

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.
# OUTPUT ERROR MISO SYSTEM IDENTIFICATION
USING FRACTIONAL MODELS

Abir Mayoufi¹,², Stéphane Victor¹, Manel Chetoui², Rachid Malti¹, Mohamed Aoun²

## Abstract

This paper deals with system identification for continuous-time multiple-input single-output (MISO) fractional differentiation models. An output error optimization algorithm is proposed for estimating all parameters, namely the coefficients and the differentiation orders. Given the high number of parameters to be estimated, the output error method can converge to a local minimum. Therefore, an initialization procedure is proposed to help the convergence to the optimum by using three variants of the algorithm. Moreover, a new definition of structured-commensurability (or S-commensurability) has been introduced to cope with the differentiation order estimation. First, a global S-commensurate order is estimated for all subsystems. Then, local S-commensurate orders are estimated (one for each subsystem). Finally the S-commensurability constraint being released, all differentiation orders are further adjusted. Estimating a global S-commensurate order greatly reduces the number of parameters and helps initializing the second variant, where local S-commensurate orders are estimated which, in turn, are used as a good initial hit for the last variant. It is known that such an initialization procedure progressively increases the number of parameters and provides good efficiency of the optimization algorithm. Monte Carlo simulation analysis are provided to evaluate the performances of this algorithm.

*MSC classifications:* 93B30, 26A33, 90C31, 49M15, 90C53, 35R11, 34A08, 37B35

*Key Words and Phrases:* system identification, continuous-time, output error, multiple-inputs and single-output, MISO system, order optimization, fractional order model.

## 1. Introduction

Fractional systems are widely used in many applicative fields for modeling diffusive phenomena since the 1950s: thermal engineering [23], electrochemistry [15, 40], biological [19, 38], electrical coil modeling [1], path
OUTPUT ERROR MISO SYSTEM IDENTIFICATION
USING FRACTIONAL MODELS

Abir Mayoufi¹,², Stéphane Victor¹, Manel Chetoui², Rachid Malti¹, Mohamed Aoun²

**Abstract**

This paper deals with system identification for continuous-time multiple-input single-output (MISO) fractional differentiation models. An output error optimization algorithm is proposed for estimating all parameters, namely the coefficients and the differentiation orders. Given the high number of parameters to be estimated, the output error method can converge to a local minimum. Therefore, an initialization procedure is proposed to help the convergence to the optimum by using three variants of the algorithm. Moreover, a new definition of structured-commensurability (or S-commensurability) has been introduced to cope with the differentiation order estimation. First, a global S-commensurate order is estimated for all subsystems. Then, local S-commensurate orders are estimated (one for each subsystem). Finally the S-commensurability constraint being released, all differentiation orders are further adjusted. Estimating a global S-commensurate order greatly reduces the number of parameters and helps initializing the second variant, where local S-commensurate orders are estimated which, in turn, are used as a good initial hit for the last variant. It is known that such an initialization procedure progressively increases the number of parameters and provides good efficiency of the optimization algorithm. Monte Carlo simulation analysis are provided to evaluate the performances of this algorithm.

*MSC classifications:* 93B30, 26A33, 90C31, 49M15, 90C53, 35R11, 34A08, 37B35

*Key Words and Phrases:* system identification, continuous-time, output error, multiple-inputs and single-output, MISO system, order optimization, fractional order model.

# 1. Introduction

Fractional systems are widely used in many applicative fields for modeling diffusive phenomena since the 1950s: thermal engineering [23], electrochemistry [15, 40], biological [19, 38], electrical coil modeling [1], path
planning through flatness principles [36, 37] path planning with fractional potential fields [28, 6] etc. Thanks to its global characterization, the fractional (or non-integer) differentiation operator is well-suited for modeling diffusion processes as opposed to the classic integer derivatives.

From a set of input and output data, a mathematical model can be obtained to describe the input/output behavior of a dynamical system by using system identification. First works on system identification were initially developed in discrete time (DT) ([12, 42, 18, 31, 39]), such as the output error (OE) method [4, 30, 29]. As dynamical systems are generally described by differential equations, system identification has been later expressed in continuous time (CT) in the preceding years [10], [43], [25], [11].

System identification with fractional order models has been widely studied: subspace method for pseudo-state-space identification [33], multi-variable non commensurate fractional systems [20], fractional order model (FOM) with time delays [22], errors-in-variables context [5], modulating functions [41, 3], inversion mechanism of functional extrema model via the differential evolution algorithms ([8]), using block pulse functions [32], measurement noise compensation based on polynomial modulating function [9], recursive identification method for fractional time-delay systems using a DT model [14], FOM identification using enhanced response sensitivity approach [17], fractional order estimation of subdiffusion equations [2].

In equation error methods [7], the estimated output is obtained from a prediction model which depends linearly on the system parameters and the delayed observations of the input and the output signals. On the contrary, the output error method is based on minimizing a criterion defined by the quadratic output error, namely the error signal between the true system output and the estimated one, by assuming that both systems are excited with the same input. Consequently, the estimated output is nonlinear in its parameters. A gradient based algorithm is used to minimize the criterion. Many fractional system identification problems have already been solved by using the output error method [24, 35, 27, 16]. All these papers concern fractional single-input single-output (SISO) system identification and no work deals with the fractional multiple-input single output (MISO) case. Also, estimating all the parameters of a fractional system turns out to estimate both coefficients and differentiation orders. Estimating the differentiation orders is tricky as the model order changes and the optimization problem is nonlinear when minimizing the output error criterion. The estimation is even more delicate when working with fractional MISO systems as the total number of parameters grows high. Hence, the main contribution of the paper is to propose an original method for initializing the gradient based algorithm by considering first a global S-commensurate
planning through flatness principles [36, 37] path planning with fractional potential fields [28, 6] etc. Thanks to its global characterization, the fractional (or non-integer) differentiation operator is well-suited for modeling diffusion processes as opposed to the classic integer derivatives.

From a set of input and output data, a mathematical model can be obtained to describe the input/output behavior of a dynamical system by using system identification. First works on system identification were initially developed in discrete time (DT) ([12, 42, 18, 31, 39]), such as the output error (OE) method [4, 30, 29]. As dynamical systems are generally described by differential equations, system identification has been later expressed in continuous time (CT) in the preceding years [10], [43], [25], [11].

System identification with fractional order models has been widely studied: subspace method for pseudo-state-space identification [33], multi-variable non commensurate fractional systems [20], fractional order model (FOM) with time delays [22], errors-in-variables context [5], modulating functions [41, 3], inversion mechanism of functional extrema model via the differential evolution algorithms ([8]), using block pulse functions [32], measurement noise compensation based on polynomial modulating function [9], recursive identification method for fractional time-delay systems using a DT model [14], FOM identification using enhanced response sensitivity approach [17], fractional order estimation of subdiffusion equations [2].

In equation error methods [7], the estimated output is obtained from a prediction model which depends linearly on the system parameters and the delayed observations of the input and the output signals. On the contrary, the output error method is based on minimizing a criterion defined by the quadratic output error, namely the error signal between the true system output and the estimated one, by assuming that both systems are excited with the same input. Consequently, the estimated output is nonlinear in its parameters. A gradient based algorithm is used to minimize the criterion. Many fractional system identification problems have already been solved by using the output error method [24, 35, 27, 16]. All these papers concern fractional single-input single-output (SISO) system identification and no work deals with the fractional multiple-input single output (MISO) case. Also, estimating all the parameters of a fractional system turns out to estimate both coefficients and differentiation orders. Estimating the differentiation orders is tricky as the model order changes and the optimization problem is nonlinear when minimizing the output error criterion. The estimation is even more delicate when working with fractional MISO systems as the total number of parameters grows high. Hence, the main contribution of the paper is to propose an original method for initializing the gradient based algorithm by considering first a global S-commensurate
order, then local S-commensurate orders, and finally by releasing the S-commensurability constraint.

The paper is organized as follows. Beginning with the problem formulation in Section 2, Section 3 presents the OE algorithm extended to MISO systems for the parameter vector estimation. Then an example illustrates the efficiency of the proposed algorithm with the three variants in Section 4. Finally, a conclusion is drawn in Section 5.

## 2. Problem formulation

2.1. **Fractional calculus.** The differentiation to an arbitrary order $\nu \in \mathbb{R}^+$ of a function $f(t)$, in the sense of Grünwald-Letnikov, is defined by:

$$p^{\nu} f(t) = \lim_{l \to 0} \frac{1}{l^{\nu}} \sum_{h=0}^{\lfloor \frac{t}{l} \rfloor} (-1)^h \binom{\nu}{h} f(t - hl) \quad (2.1)$$

where $\lfloor \cdot \rfloor$ is the floor operator and $\binom{\nu}{h}$ is the Newton's binomial generalized to fractional orders as follows:

$$\binom{\nu}{h} = \frac{\Gamma(\nu + 1)}{\Gamma(h + 1)\Gamma(\nu - h + 1)} = \frac{\nu(\nu - 1)\dots(\nu - h + 1)}{h!}.$$

Consequently, the fractional derivative of a function has a global characterization as the whole past of the function is required.

For $\nu = 1$, all binomials $\binom{\nu}{h} = 0$ when $h - \nu \in \mathbb{N}$, thus leading equation (2.1) to the classic definition of integer order derivative:

$$pf(t) = \lim_{l \to 0} \frac{f(t) - f(t - l)}{l}. \quad (2.2)$$

For numerical evaluation of the $\nu$-th fractional derivative, parameter $l$ in (2.1) is replaced by a sampling period $T$ and the limit is dropped:

$$p^{\nu} f(t) = \frac{1}{T^{\nu}} \sum_{h=0}^{\lfloor \frac{t}{T} \rfloor} (-1)^h \binom{\nu}{h} f(t - Th) + O(T). \quad (2.3)$$

So defined, as stated in [23, section 7.4], the error term $O(T)$ is proportional to the sampling period. Therefore, to make the approximation error negligible, the sampling period must be sufficiently small.

<footnote>Note that a too small sampling time may also provide numerical problems such as stability in digital implementations. In this case, discrete rational approximations may be used.</footnote>
order, then local S-commensurate orders, and finally by releasing the S-commensurability constraint.

The paper is organized as follows. Beginning with the problem formulation in Section 2, Section 3 presents the OE algorithm extended to MISO systems for the parameter vector estimation. Then an example illustrates the efficiency of the proposed algorithm with the three variants in Section 4. Finally, a conclusion is drawn in Section 5.

## 2. Problem formulation

2.1. **Fractional calculus.** The differentiation to an arbitrary order $\nu \in \mathbb{R}^+$ of a function $f(t)$, in the sense of Grünwald-Letnikov, is defined by:

$$p^{\nu} f(t) = \lim_{l \to 0} \frac{1}{l^{\nu}} \sum_{h=0}^{\lfloor \frac{t}{l} \rfloor} (-1)^h \binom{\nu}{h} f(t - hl) \quad (2.1)$$

where $\lfloor \cdot \rfloor$ is the floor operator and $\binom{\nu}{h}$ is the Newton's binomial generalized to fractional orders as follows:

$$\binom{\nu}{h} = \frac{\Gamma(\nu + 1)}{\Gamma(h + 1)\Gamma(\nu - h + 1)} = \frac{\nu(\nu - 1)\dots(\nu - h + 1)}{h!}.$$

Consequently, the fractional derivative of a function has a global characterization as the whole past of the function is required.

For $\nu = 1$, all binomials $\binom{\nu}{h} = 0$ when $h - \nu \in \mathbb{N}$, thus leading equation (2.1) to the classic definition of integer order derivative:

$$pf(t) = \lim_{l \to 0} \frac{f(t) - f(t - l)}{l}. \quad (2.2)$$

For numerical evaluation of the $\nu$-th fractional derivative, parameter $l$ in (2.1) is replaced by a sampling period $T$ and the limit is dropped:

$$p^{\nu} f(t) = \frac{1}{T^{\nu}} \sum_{h=0}^{\lfloor \frac{t}{T} \rfloor} (-1)^h \binom{\nu}{h} f(t - Th) + O(T). \quad (2.3)$$

So defined, as stated in [23, section 7.4], the error term $O(T)$ is proportional to the sampling period. Therefore, to make the approximation error negligible, the sampling period must be sufficiently small.

<footnote>Note that a too small sampling time may also provide numerical problems such as stability in digital implementations. In this case, discrete rational approximations may be used.</footnote>
The Laplace transform of the $\nu$-th derivative of a causal function $f$ ($f(t) = 0 \forall t \le 0$), is given by [23, section 2.8.4]:

$$ \mathcal{L}\{p^{\nu}x(t)\} = s^{\nu}X(s), \qquad (2.4) $$

$s$ denoting the Laplace variable. This result is in accordance with the classic rational case when $\nu$ is an integer.

## 2.2. MISO fractional order system.
A MISO fractional order system (FOS) such as represented in Fig.1, can be described as the following:

$$ S: \left\{ \begin{array}{l} y_k(t) = G_k(p)u_k(t) \\ y(t) = \sum_{k=1}^{K} y_k(t) \\ y^*(t_h) = y(t_h) + e(t_h) \end{array} \right. \qquad (2.5) $$

where $K \in \mathbb{N}^*$ represents the number of inputs, and consequently the number of subsystems. The vector $u(t) = [u_1(t), \dots, u_K(t)]$ contains the uncorrelated input signals and $p$ denotes the differential operator $p = \frac{d}{dt}$. The input-output data $(u(t), y(t))$ are collected at regular samples and assumed large enough to guarantee convergence of the estimated parameters to the true ones. The quasi-stationary input signals $\{u_k(t), 0 \le t \le t_{final}, k = 1, \dots, K\}$ applied to the MISO system are assumed persistently exciting and give rise to the output signal $\{y(t), 0 \le t \le t_{final}\}$. The output measurement $y(t)$ is assumed corrupted by a discrete-time white noise $e(t_h)$ and $G_k$ is the fractional operator that relates the input signal $u_k(t)$ to its noise-free output $y_k(t)$:

$$ G_k(p) = \frac{B_k(p)}{A_k(p)} = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{\beta_{j,k}}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{\alpha_{i,k}}} \qquad (2.6) $$

where $(a_{i,k}, b_{j,k}) \in \mathbb{R}^2$, for all $i, j$ and $k$ and the differentiation orders $(\alpha_{i,k}, \beta_{j,k})$ are ordered positive real numbers for identifiability purpose:

$$ \left\{ \begin{array}{l} 0 < \alpha_{1,k} < \alpha_{2,k} < \dots < \alpha_{N_k,k} \quad \forall k = 1, \dots, K. \\ 0 < \beta_{0,k} < \beta_{1,k} < \dots < \beta_{M_k,k} \end{array} \right. $$

Let be strictly proper operators $G_k(p)$ with $\alpha_{N_k,k} > \beta_{M_k,k}, \forall k$. Also, $B_k$ and $A_k$ are coprime polynomials, and all the transfer functions $G_k$ are asymptotically stable [21, 26].
The Laplace transform of the $\nu$-th derivative of a causal function $f$ ($f(t) = 0 \forall t \le 0$), is given by [23, section 2.8.4]:

$$ \mathcal{L}\{p^{\nu}x(t)\} = s^{\nu}X(s), \qquad (2.4) $$

$s$ denoting the Laplace variable. This result is in accordance with the classic rational case when $\nu$ is an integer.

**2.2. MISO fractional order system.** A MISO fractional order system (FOS) such as represented in Fig.1, can be described as the following:

$$ S: \left\{ \begin{array}{l} y_k(t) = G_k(p)u_k(t) \\ y(t) = \sum_{k=1}^{K} y_k(t) \\ y^*(t_h) = y(t_h) + e(t_h) \end{array} \right. \qquad (2.5) $$

where $K \in \mathbb{N}^*$ represents the number of inputs, and consequently the number of subsystems. The vector $u(t) = [u_1(t), \dots, u_K(t)]$ contains the uncorrelated input signals and $p$ denotes the differential operator $p = \frac{d}{dt}$. The input-output data $(u(t), y(t))$ are collected at regular samples and assumed large enough to guarantee convergence of the estimated parameters to the true ones. The quasi-stationary input signals $\{u_k(t), 0 \le t \le t_{final}, k = 1, \dots, K\}$ applied to the MISO system are assumed persistently exciting and give rise to the output signal $\{y(t), 0 \le t \le t_{final}\}$. The output measurement $y(t)$ is assumed corrupted by a discrete-time white noise $e(t_h)$ and $G_k$ is the fractional operator that relates the input signal $u_k(t)$ to its noise-free output $y_k(t)$:

$$ G_k(p) = \frac{B_k(p)}{A_k(p)} = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{\beta_{j,k}}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{\alpha_{i,k}}} \qquad (2.6) $$

where $(a_{i,k}, b_{j,k}) \in \mathbb{R}^2$, for all $i, j$ and $k$ and the differentiation orders $(\alpha_{i,k}, \beta_{j,k})$ are ordered positive real numbers for identifiability purpose:

$$ \left\{ \begin{array}{l} 0 < \alpha_{1,k} < \alpha_{2,k} < \dots < \alpha_{N_k,k} \quad \forall k = 1, \dots, K. \\ 0 < \beta_{0,k} < \beta_{1,k} < \dots < \beta_{M_k,k} \end{array} \right. $$

Let be strictly proper operators $G_k(p)$ with $\alpha_{N_k,k} > \beta_{M_k,k}, \forall k$. Also, $B_k$ and $A_k$ are coprime polynomials, and all the transfer functions $G_k$ are asymptotically stable [21, 26].
<figure><img src="image_1.png" /><figcaption>FIGURE 1. Fractional MISO system structure</figcaption></figure>

DEFINITION 2.1 (S-commensurability). A SISO system $G_1$, for a given model structure (2.6), characterized by the number of terms $M_1$ and $N_1$, has a structured commensurate (or S-commensurate) order $\nu = \alpha_{1,1}$, if all its differentiation orders are successive integer multiples of $\nu$:

$$G_1(p) = \frac{\sum_{j=0}^{M_1} b_{j,1} p^{j\nu}}{1 + \sum_{i=1}^{N_1} a_{i,1} p^{i\nu}}. \qquad (2.7)$$

For example, the following transfer function:

$$G_1(s) = \frac{1}{1 + 1.5s^{1.2} + 0.5s^{2.8}}$$

is commensurate of order 0.4 according to [21], namely:

$$G_1(s) = \frac{1}{1 + 1.5s^{3 \times 0.4} + 0.5s^{6 \times 0.4}}$$

According to the new definition 2.1 of S-commensurability, $G_1$ is not S-commensurate, because the differentiation orders are not successive integer multiples of $\alpha_{1,1} = 1.2$.

Note that the number of terms in the numerator and the denominator $M_1$ and $N_1$ do not change from (2.6) to (2.7).

DEFINITION 2.2 (Local S-commensurability). A MISO system, for a given model structure (2.6), has local S-commensurate orders $\nu_k$ if each

<footnote>In the classic case ([21]), the definition of a commensurate system is independent of the model structure and may generate very high number of parameters when the commensurate order is very low. The S-commensurability definition is linked to a model structure and therefore fixes the number of parameters to a prescribed value.</footnote>
<figure><img src="image_1.png" /><figcaption>FIGURE 1. Fractional MISO system structure</figcaption></figure>

DEFINITION 2.1 (S-commensurability). A SISO system $G_1$, for a given model structure (2.6), characterized by the number of terms $M_1$ and $N_1$, has a structured commensurate (or S-commensurate) order $\nu = \alpha_{1,1}$, if all its differentiation orders are successive integer multiples of $\nu$:

$$G_1(p) = \frac{\sum_{j=0}^{M_1} b_{j,1} p^{j\nu}}{1 + \sum_{i=1}^{N_1} a_{i,1} p^{i\nu}}. \qquad (2.7)$$

For example, the following transfer function:

$$G_1(s) = \frac{1}{1 + 1.5s^{1.2} + 0.5s^{2.8}}$$

is commensurate of order 0.4 according to [21], namely:

$$G_1(s) = \frac{1}{1 + 1.5s^{3 \times 0.4} + 0.5s^{6 \times 0.4}}$$

According to the new definition 2.1 of S-commensurability, $G_1$ is not S-commensurate, because the differentiation orders are not successive integer multiples of $\alpha_{1,1} = 1.2$.

Note that the number of terms in the numerator and the denominator $M_1$ and $N_1$ do not change from (2.6) to (2.7).

DEFINITION 2.2 (Local S-commensurability). A MISO system, for a given model structure (2.6), has local S-commensurate orders $\nu_k$ if each

<footnote>In the classic case ([21]), the definition of a commensurate system is independent of the model structure and may generate very high number of parameters when the commensurate order is very low. The S-commensurability definition is linked to a model structure and therefore fixes the number of parameters to a prescribed value.</footnote>
subsystem $G_k$ has its own S-commensurate order (see Definition 2.1)

$$ \nu_k = \alpha_{1,k} \quad \text{for} \quad k = 1, \dots, K. $$

Therefore, each transfer function can be rewritten as:

$$ G_k(p) = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{j\nu_k}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{i\nu_k}}. \qquad (2.8) $$

**DEFINITION 2.3 (Global S-commensurability).** A MISO system, for a given model structure (2.6), has a global S-commensurate order $\nu$ if all its subsystems have the same S-commensurate order according to Definition 2.1:

$$ \nu = \alpha_{1,1} = \alpha_{1,2} = \dots = \alpha_{1,K}, $$

Therefore, each transfer function can be rewritten as:

$$ G_k(p) = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{j\nu}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{i\nu}}. \qquad (2.9) $$

Again, note that the numerator and denominator orders $M_k$ and $N_k$ do not change from (2.6) to (2.9).

For a S-commensurate FOS, the most used stability theorem is the one proposed by [21] extended to commensurate orders between 0 and 2.

**THEOREM 2.1.** *Stability theorem*

Let $G_1$ be a SISO S-commensurate transfer function and $\nu$ its S-commensurate order. $G_1(s) = \frac{Q_\nu(s)}{P_\nu(s)}$ is stable in the bounded input bounded output sense if and only if:

$$ 0 < \nu < 2 \qquad (2.10) $$

and

$$ \forall s_k \in \mathbb{C}, P_\nu(s_k) = 0 \quad \text{such as} \quad |\arg(s_k)| > \nu \frac{\pi}{2}. \qquad (2.11) $$

□

If all subsystems are stable then the MISO fractional system is stable too.
subsystem $G_k$ has its own S-commensurate order (see Definition 2.1)

$$ \nu_k = \alpha_{1,k} \quad \text{for} \quad k = 1, \dots, K. $$

Therefore, each transfer function can be rewritten as:

$$ G_k(p) = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{j\nu_k}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{i\nu_k}}. \qquad (2.8) $$

**DEFINITION 2.3 (Global S-commensurability).** A MISO system, for a given model structure (2.6), has a global S-commensurate order $\nu$ if all its subsystems have the same S-commensurate order according to Definition 2.1:

$$ \nu = \alpha_{1,1} = \alpha_{1,2} = \dots = \alpha_{1,K}, $$

Therefore, each transfer function can be rewritten as:

$$ G_k(p) = \frac{\sum_{j=0}^{M_k} b_{j,k} p^{j\nu}}{1 + \sum_{i=1}^{N_k} a_{i,k} p^{i\nu}}. \qquad (2.9) $$

Again, note that the numerator and denominator orders $M_k$ and $N_k$ do not change from (2.6) to (2.9).

For a S-commensurate FOS, the most used stability theorem is the one proposed by [21] extended to commensurate orders between 0 and 2.

**THEOREM 2.1.** *Stability theorem*

Let $G_1$ be a SISO S-commensurate transfer function and $\nu$ its S-commensurate order. $G_1(s) = \frac{Q_\nu(s)}{P_\nu(s)}$ is stable in the bounded input bounded output sense if and only if:

$$ 0 < \nu < 2 \qquad (2.10) $$

and

$$ \forall s_k \in \mathbb{C}, P_\nu(s_k) = 0 \quad \text{such as} \quad |\arg(s_k)| > \nu \frac{\pi}{2}. \qquad (2.11) $$

□

If all subsystems are stable then the MISO fractional system is stable too.
2.3. **Problem formulation.** The aim is to estimate all the parameters of fractional MISO model $\mathcal{S}$ described by equation (2.5) by using $H$ samples of input/output data $\{u_1(t_h), \dots, u_K(t_h), y^*(t_h)\}_{h=1}^H$. Note that in conformity with section 2.2, $t_1 = 0$ and $t_H = t_{final}$.

The parameter vector $\theta$ of the MISO system (2.5) is defined as

$$ \theta = \begin{bmatrix} \rho \\ \mu \end{bmatrix} \qquad (2.12) $$

• where $\rho$ has all the coefficients

$$ \rho = [\rho_1, \dots, \rho_K]^T \qquad (2.13) $$

with

$$ \rho_k = [b_{0,k}, b_{1,k}, \dots, b_{M_k,k}, a_{1,k}, \dots, a_{N_k,k}] $$

for $k = 1, \dots, K$; hence, $\rho$ has a total number of coefficients of $\sum_{k=1}^K (N_k + M_k + 1)$

• and where $\mu$ has all the differentiation orders, that can be defined according to the following cases.

Case 1: if a global S-commensurate order $\nu$ is sought, then the differentiation order vector is reduced to

$$ \mu = \nu; \qquad (2.14) $$

Case 2: if local S-commensurate orders $\nu_k$ are sought, then $\mu$ is extended to $K$ parameters

$$ \mu = [\nu_1, \dots, \nu_K]^T; \qquad (2.15) $$

Case 3: if the MISO model is non commensurate, then $\mu$ gathers all the differentiation orders

$$ \mu = [\mu_1, \dots, \mu_K]^T \qquad (2.16) $$

with

$$ \mu_k = [\beta_{0,k}, \beta_{1,k}, \dots, \beta_{M_k,k}, \alpha_{1,k}, \dots, \alpha_{N_k,k}] $$

for $k = 1, \dots, K$; hence, $\mu$ has a total number of differentiation orders of $\sum_{k=1}^K (N_k + M_k + 1)$.

All the proposed system identification algorithms are meant to be independent from time-domain simulation algorithms. Nonetheless, FOS should be correctly simulated with negligible simulation errors so that the parameter estimation is consistent.
2.3. **Problem formulation.** The aim is to estimate all the parameters of fractional MISO model $\mathcal{S}$ described by equation (2.5) by using $H$ samples of input/output data $\{u_1(t_h), ..., u_K(t_h), y^*(t_h)\}_{h=1}^H$. Note that in conformity with section 2.2, $t_1 = 0$ and $t_H = t_{final}$.

The parameter vector $\theta$ of the MISO system (2.5) is defined as

$$ \theta = \begin{bmatrix} \rho \\ \mu \end{bmatrix} \qquad (2.12) $$

• where $\rho$ has all the coefficients

$$ \rho = [\rho_1, \dots, \rho_K]^T \qquad (2.13) $$

with

$$ \rho_k = [b_{0,k}, b_{1,k}, \dots, b_{M_k,k}, a_{1,k}, \dots, a_{N_k,k}] $$

for $k = 1, \dots, K$; hence, $\rho$ has a total number of coefficients of $\sum_{k=1}^K (N_k + M_k + 1)$

• and where $\mu$ has all the differentiation orders, that can be defined according to the following cases.

Case 1: if a global S-commensurate order $\nu$ is sought, then the differentiation order vector is reduced to

$$ \mu = \nu; \qquad (2.14) $$

Case 2: if local S-commensurate orders $\nu_k$ are sought, then $\mu$ is extended to $K$ parameters

$$ \mu = [\nu_1, \dots, \nu_K]^T; \qquad (2.15) $$

Case 3: if the MISO model is non commensurate, then $\mu$ gathers all the differentiation orders

$$ \mu = [\mu_1, \dots, \mu_K]^T \qquad (2.16) $$

with

$$ \mu_k = [\beta_{0,k}, \beta_{1,k}, \dots, \beta_{M_k,k}, \alpha_{1,k}, \dots, \alpha_{N_k,k}] $$

for $k = 1, \dots, K$; hence, $\mu$ has a total number of differentiation orders of $\sum_{k=1}^K (N_k + M_k + 1)$.

All the proposed system identification algorithms are meant to be independent from time-domain simulation algorithms. Nonetheless, FOS should be correctly simulated with negligible simulation errors so that the parameter estimation is consistent.
### 3. Output Error method for CT fractional MISO models

The true model structure being assumed known, the parameter vector $\theta$ is estimated using a nonlinear optimization method, based on Levenberg-Marquardt algorithm.

The estimation problem is formulated as a minimization problem of the $l_2$-norm:

$$J(\hat{\theta}) = \frac{1}{2} \|\varepsilon(t, \hat{\theta})\|_2, \qquad (3.1)$$

where the output error $\varepsilon(t, \hat{\theta})$ is defined by:

$$\varepsilon(t, \hat{\theta}) = y^*(t) - \hat{y}(t) \qquad (3.2)$$

and the estimated output $\hat{y}(t)$ is defined by:

$$\hat{y}(t) = \sum_{k=1}^{K} \hat{y}_k(t) \qquad (3.3)$$

with $\hat{y}_k(t) = \hat{G}_k(p)u_k(t)$ and $\hat{G}_k(p)$ denote the estimated model.

The algorithm being iterative, therefore the parameter vector $\theta$ in (2.12) is denoted by $\theta^{iter}$. The OE algorithm for fractional MISO system is described in Algorithm 1.

After convergence, the parametric covariance matrix can be computed from the approximated Hessian $\mathcal{H}$, namely:

$$\mathcal{H} = \frac{\partial \varepsilon^T}{\partial \theta} \frac{\partial \varepsilon}{\partial \theta} \qquad (3.8)$$

where $\frac{\partial \varepsilon}{\partial \theta}$ is the error sensitivity function given by:

$$\frac{\partial \varepsilon}{\partial \theta} = \frac{\partial \varepsilon}{\partial [\rho^T \mu^T]^T}. \qquad (3.9)$$

The coefficient error sensitivity is given by:

$$\frac{\partial \varepsilon}{\partial \rho} = \left[ \frac{\partial \varepsilon}{\partial \rho_1}, \dots, \frac{\partial \varepsilon}{\partial \rho_K} \right], \qquad (3.10)$$

where

$$\frac{\partial \varepsilon}{\partial \rho_k} = -\frac{\partial \hat{y}_k}{\partial \rho_k} = -\left[ \frac{\partial \hat{y}_k}{\partial b_{0,k}}, \dots, \frac{\partial \hat{y}_k}{\partial b_{M_k,k}}, \frac{\partial \hat{y}_{u_k}}{\partial a_{1,k}}, \dots, \frac{\partial \hat{y}_{u_k}}{\partial a_{N_k,k}} \right] \qquad (3.11)$$
### 3. Output Error method for CT fractional MISO models

The true model structure being assumed known, the parameter vector $\theta$ is estimated using a nonlinear optimization method, based on Levenberg-Marquardt algorithm.

The estimation problem is formulated as a minimization problem of the $l_2$-norm:

$$J(\hat{\theta}) = \frac{1}{2} \|\varepsilon(t, \hat{\theta})\|_2, \qquad (3.1)$$

where the output error $\varepsilon(t, \hat{\theta})$ is defined by:

$$\varepsilon(t, \hat{\theta}) = y^*(t) - \hat{y}(t) \qquad (3.2)$$

and the estimated output $\hat{y}(t)$ is defined by:

$$\hat{y}(t) = \sum_{k=1}^{K} \hat{y}_k(t) \qquad (3.3)$$

with $\hat{y}_k(t) = \hat{G}_k(p)u_k(t)$ and $\hat{G}_k(p)$ denote the estimated model.

The algorithm being iterative, therefore the parameter vector $\theta$ in (2.12) is denoted by $\theta^{iter}$. The OE algorithm for fractional MISO system is described in Algorithm 1.

After convergence, the parametric covariance matrix can be computed from the approximated Hessian $\mathcal{H}$, namely:

$$\mathcal{H} = \frac{\partial \varepsilon^T}{\partial \theta} \frac{\partial \varepsilon}{\partial \theta} \qquad (3.8)$$

where $\frac{\partial \varepsilon}{\partial \theta}$ is the error sensitivity function given by:

$$\frac{\partial \varepsilon}{\partial \theta} = \frac{\partial \varepsilon}{\partial [\rho^T \mu^T]^T} \qquad (3.9)$$

The coefficient error sensitivity is given by:

$$\frac{\partial \varepsilon}{\partial \rho} = \left[ \frac{\partial \varepsilon}{\partial \rho_1}, \dots, \frac{\partial \varepsilon}{\partial \rho_K} \right], \qquad (3.10)$$

where

$$\frac{\partial \varepsilon}{\partial \rho_k} = -\frac{\partial \hat{y}_k}{\partial \rho_k} = -\left[ \frac{\partial \hat{y}_k}{\partial b_{0,k}}, \dots, \frac{\partial \hat{y}_k}{\partial b_{M_k,k}}, \frac{\partial \hat{y}_{u_k}}{\partial a_{1,k}}, \dots, \frac{\partial \hat{y}_{u_k}}{\partial a_{N_k,k}} \right] \qquad (3.11)$$
<caption>**Algorithm 1: MISO OE algorithm**</caption>

**Step 1: Initialization**
iter = 0
Initialize the parameter vector $\theta^0 = [\rho^0, \mu^0]^T$
Compute the cost function $J(\theta^0)$ from (3.1)

**Step 2: Iterative Levenberg-Marquardt algorithm**
do

i. iter = iter + 1

ii. Parameter estimation
Compute the parameter vector $\hat{\theta}^{\text{iter}}$ as

$$ \hat{\theta}^{\text{iter}+1} = \hat{\theta}^{\text{iter}} - \left\{ [\mathcal{H} + \xi I]^{-1} \frac{\partial J}{\partial \hat{\theta}} \right\|_{\hat{\theta}^{\text{iter}}}, \quad (3.4) $$

with

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial J}{\partial \hat{\theta}} = \sum_{h=1}^{H} \frac{\partial \varepsilon(t_h)^T}{\partial \hat{\theta}} \varepsilon(t_h) : \text{gradient} \\[1.5em] \displaystyle \mathcal{H} \approx \sum_{h=1}^{H} \frac{\partial \varepsilon(t_h)^T}{\partial \hat{\theta}} \frac{\partial \varepsilon(t_h)}{\partial \hat{\theta}} : \text{pseudo-Hessian} \\[1.5em] \displaystyle \xi : \text{Marquardt parameter} \end{array} \right. \quad (3.5) $$

while $\left| \frac{\hat{\theta}_{\ell}^{\text{iter}} - \hat{\theta}_{\ell}^{\text{iter}-1}}{\hat{\theta}_{\ell}^{\text{iter}-1}} \right| > \epsilon$

where $\hat{\theta}_{\ell}^{\text{iter}}$ corresponds to the $\ell$-th element of the vector $\hat{\theta}^{\text{iter}}$ at iteration iter.

**Step3: Covariance matrix estimation**

Compute the parametric covariance matrix:

$$ \hat{P} = \sigma_{opt}^2 \mathcal{H}^{-1} \quad (3.6) $$

where $\mathcal{H}$ denotes the approximated Hessian defined in (3.8) and $\hat{\sigma}^2$ is the variance of the empirical estimation of the noise, $\hat{e}(t)$, given by:

$$ \hat{e}(t) = y^*(t) - \sum_{k=1}^{K} y_k^{\text{iter}}(t) \quad (3.7) $$

with $y_k^{\text{iter}}(t)$ being calculated at the last iteration.

for $k = 1, ..., K$. The output sensitivity functions $\frac{\partial \hat{y}_k}{\partial b_{j,k}}$ and $\frac{\partial \hat{y}_k}{\partial a_{i,k}}$ are given by:

$$ \frac{\partial \hat{y}_k}{\partial b_{j,k}} = \frac{p^{\hat{\beta}_j}}{1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{\hat{\alpha}_i}} u_k(t), \quad \forall j = 0, ..., M_k \quad (3.12) $$
<caption>**Algorithm 1: MISO OE algorithm**</caption>

**Step 1: Initialization**
iter = 0
Initialize the parameter vector $\theta^0 = [\rho^0, \mu^0]^T$
Compute the cost function $J(\theta^0)$ from (3.1)

**Step 2: Iterative Levenberg-Marquardt algorithm**
do

i. iter = iter + 1

ii. Parameter estimation
Compute the parameter vector $\hat{\theta}^{\text{iter}}$ as

$$ \hat{\theta}^{\text{iter}+1} = \hat{\theta}^{\text{iter}} - \left\{ [\mathcal{H} + \xi I]^{-1} \frac{\partial J}{\partial \hat{\theta}} \right\|_{\hat{\theta}^{\text{iter}}}, \quad (3.4) $$

with

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial J}{\partial \hat{\theta}} = \sum_{h=1}^{H} \frac{\partial \varepsilon(t_h)^T}{\partial \hat{\theta}} \varepsilon(t_h) : \text{gradient} \\[1em] \displaystyle \mathcal{H} \approx \sum_{h=1}^{H} \frac{\partial \varepsilon(t_h)^T}{\partial \hat{\theta}} \frac{\partial \varepsilon(t_h)}{\partial \hat{\theta}} : \text{pseudo-Hessian} \\[1em] \displaystyle \xi : \text{Marquardt parameter} \end{array} \right. \quad (3.5) $$

while $\left| \frac{\hat{\theta}_{\ell}^{\text{iter}} - \hat{\theta}_{\ell}^{\text{iter}-1}}{\hat{\theta}_{\ell}^{\text{iter}-1}} \right| > \epsilon$
where $\hat{\theta}_{\ell}^{\text{iter}}$ corresponds to the $\ell$-th element of the vector $\hat{\theta}^{\text{iter}}$ at iteration iter.

**Step3: Covariance matrix estimation**
Compute the parametric covariance matrix:

$$ \hat{P} = \sigma_{opt}^2 \mathcal{H}^{-1} \quad (3.6) $$

where $\mathcal{H}$ denotes the approximated Hessian defined in (3.8) and $\hat{\sigma}^2$ is the variance of the empirical estimation of the noise, $\hat{e}(t)$, given by:

$$ \hat{e}(t) = y^*(t) - \sum_{k=1}^{K} y_k^{\text{iter}}(t) \quad (3.7) $$

with $y_k^{\text{iter}}(t)$ being calculated at the last iteration.

for $k = 1, ..., K$. The output sensitivity functions $\frac{\partial \hat{y}_k}{\partial b_{j,k}}$ and $\frac{\partial \hat{y}_k}{\partial a_{i,k}}$ are given by:

$$ \frac{\partial \hat{y}_k}{\partial b_{j,k}} = \frac{p^{\hat{\beta}_j}}{1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{\hat{\alpha}_i}} u_k(t), \quad \forall j = 0, ..., M_k \quad (3.12) $$
and

$$ \frac{\partial \hat{y}_k}{\partial a_{i,k}} = - \frac{\sum_{j=0}^{M_k} \hat{b}_{j,k} p^{\hat{\beta}_j + \hat{\alpha}_i}}{\left(1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{\hat{\alpha}_i}\right)^2} u_k(t), \quad \forall i = 1, \dots, N_k. \quad (3.13) $$

The differentiation order error sensitivity functions are defined according to $\mu$:

Case 1: If $\mu$ is defined as a global S-commensurate order $\nu$ for all subsystems as in (2.14), then

$$ \frac{\partial \varepsilon}{\partial \mu} = \frac{\partial \varepsilon}{\partial \nu} = - \sum_{k=1}^{K} \frac{\partial \hat{y}_k}{\partial \nu} \quad (3.14) $$

where the output sensitivity function is given, for $k = 1, \dots, K$, by:

$$ \frac{\partial \hat{y}_k}{\partial \nu} = \left[ \sum_{j=0}^{M_k} j \hat{b}_{j,k} p^{j\nu} + \sum_{j=0}^{M_k} \sum_{i=1}^{N_k} (j-i) \hat{b}_{j,k} \hat{a}_{i,k} p^{(i+j)\nu} \right] \times \frac{\ln(p)}{\left(1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{i\nu}\right)^2} u_k(t). \quad (3.15) $$

Case 2: If $\mu$ is defined as the vector of local S-commensurate orders $\nu_k$, as in (2.15), then

$$ \begin{aligned} \frac{\partial \varepsilon}{\partial \mu} &= \left[ \frac{\partial \varepsilon}{\partial \nu_1}, \dots, \frac{\partial \varepsilon}{\partial \nu_K} \right] \\ &= \left[ -\frac{\partial \hat{y}_1}{\partial \nu_1}, \dots, -\frac{\partial \hat{y}_K}{\partial \nu_K} \right] \end{aligned} \quad (3.16) $$

where the output sensitivity functions are obtained by replacing $\nu$ by $\nu_k$ in (3.15), for $k = 1, \dots, K$.

Case 3: If $\mu$ is the vector of all the differentiation orders, as in (2.16), then

$$ \frac{\partial \varepsilon}{\partial \mu} = \left[ \frac{\partial \varepsilon}{\partial \mu_1}, \dots, \frac{\partial \varepsilon}{\partial \mu_K} \right], \quad (3.17) $$

with

$$ \frac{\partial \varepsilon}{\partial \mu_k} = \left[ -\frac{\partial \hat{y}_k}{\partial \beta_{0,k}}, \dots, -\frac{\partial \hat{y}_k}{\partial \beta_{M_k,k}}, -\frac{\partial \hat{y}_k}{\partial \alpha_{1,k}}, \dots, -\frac{\partial \hat{y}_k}{\partial \alpha_{N_k,k}} \right] \quad (3.18) $$
and

$$ \frac{\partial \hat{y}_k}{\partial a_{i,k}} = - \frac{\sum_{j=0}^{M_k} \hat{b}_{j,k} p^{\hat{\beta}_j + \hat{\alpha}_i}}{\left(1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{\hat{\alpha}_i}\right)^2} u_k(t), \quad \forall i = 1, \dots, N_k. \quad (3.13) $$

The differentiation order error sensitivity functions are defined according to $\mu$:

Case 1: If $\mu$ is defined as a global S-commensurate order $\nu$ for all subsystems as in (2.14), then

$$ \frac{\partial \varepsilon}{\partial \mu} = \frac{\partial \varepsilon}{\partial \nu} = - \sum_{k=1}^{K} \frac{\partial \hat{y}_k}{\partial \nu} \quad (3.14) $$

where the output sensitivity function is given, for $k = 1, \dots, K$, by:

$$ \frac{\partial \hat{y}_k}{\partial \nu} = \left[ \sum_{j=0}^{M_k} j \hat{b}_{j,k} p^{j\nu} + \sum_{j=0}^{M_k} \sum_{i=1}^{N_k} (j-i) \hat{b}_{j,k} \hat{a}_{i,k} p^{(i+j)\nu} \right] \times \frac{\ln(p)}{\left(1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{i\nu}\right)^2} u_k(t). \quad (3.15) $$

Case 2: If $\mu$ is defined as the vector of local S-commensurate orders $\nu_k$, as in (2.15), then

$$ \begin{aligned} \frac{\partial \varepsilon}{\partial \mu} &= \left[ \frac{\partial \varepsilon}{\partial \nu_1}, \dots, \frac{\partial \varepsilon}{\partial \nu_K} \right] \\ &= \left[ -\frac{\partial \hat{y}_1}{\partial \nu_1}, \dots, -\frac{\partial \hat{y}_K}{\partial \nu_K} \right] \end{aligned} \quad (3.16) $$

where the output sensitivity functions are obtained by replacing $\nu$ by $\nu_k$ in (3.15), for $k = 1, \dots, K$.

Case 3: If $\mu$ is the vector of all the differentiation orders, as in (2.16), then

$$ \frac{\partial \varepsilon}{\partial \mu} = \left[ \frac{\partial \varepsilon}{\partial \mu_1}, \dots, \frac{\partial \varepsilon}{\partial \mu_K} \right], \quad (3.17) $$

with

$$ \frac{\partial \varepsilon}{\partial \mu_k} = \left[ -\frac{\partial \hat{y}_k}{\partial \beta_{0,k}}, \dots, -\frac{\partial \hat{y}_k}{\partial \beta_{M_k,k}}, -\frac{\partial \hat{y}_k}{\partial \alpha_{1,k}}, \dots, -\frac{\partial \hat{y}_k}{\partial \alpha_{N_k,k}} \right] \quad (3.18) $$
where the output sensitivity functions are given for $k = 1, \dots, K$:

$$ \frac{\partial \hat{y}_k}{\partial \beta_{j,k}} = \ln(p) \frac{\hat{b}_{j,k} p^{\hat{\beta}_{j,k}}}{1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{\hat{\alpha}_{i,k}}} u_k(t) \quad (3.19) $$

for $j = 0, \dots, M_k$ and

$$ \frac{\partial \hat{y}_k}{\partial \alpha_{\ell,k}} = - \ln(p) \hat{a}_{\ell,k} p^{\hat{\alpha}_{\ell,k}} \frac{\sum_{j=0}^{M_k} \hat{b}_{j,k} p^{\hat{\beta}_{j,k}}}{\left(1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{\hat{\alpha}_{i,k}}\right)^2} u_k(t), \quad (3.20) $$

for $\ell = 1, \dots, N_k$.

In the output sensitivity functions (3.15), (3.19) and (3.20), the term $\ln(p)$ cannot be simulated in practice, therefore a numerical computation is used to calculate these sensitivity functions.

## 4. Simulation example

To illustrate the performances of the proposed OE method for MISO FOS as described in Algorithm 1, a numerical example is presented in this section.

Note that in real implementation, unknown initial conditions can introduce a transient of the simulated model, thus leading to bias on the estimation. Therefore, a history-function-based initialization [13] or alternatively, an infinite-dimensional state-space representation [34] can be used to better initialize fractional systems, so that the parameter estimation gets more consistent.

System identification, and particularly for fractional systems, is usually undertaken only around null initial conditions or around a working point. Therefore, the identified model is valid only for small variations around the considered working point.

To avoid local minima first, a global S-commensurate order is estimated for all subsystems which is used as an initial hit for the second variant where local S-commensurate orders are estimated, which in turn is used as an initial hit for the last variant to estimate all differentiation orders.
where the output sensitivity functions are given for $k = 1, \dots, K$:

$$ \frac{\partial \hat{y}_k}{\partial \beta_{j,k}} = \ln(p) \frac{\hat{b}_{j,k} p^{\hat{\beta}_{j,k}}}{1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{\hat{\alpha}_{i,k}}} u_k(t) \quad (3.19) $$

for $j = 0, \dots, M_k$ and

$$ \frac{\partial \hat{y}_k}{\partial \alpha_{\ell,k}} = - \ln(p) \hat{a}_{\ell,k} p^{\hat{\alpha}_{\ell,k}} \frac{\sum_{j=0}^{M_k} \hat{b}_{j,k} p^{\hat{\beta}_{j,k}}}{\left(1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} p^{\hat{\alpha}_{i,k}}\right)^2} u_k(t), \quad (3.20) $$

for $\ell = 1, \dots, N_k$.

In the output sensitivity functions (3.15), (3.19) and (3.20), the term $\ln(p)$ cannot be simulated in practice, therefore a numerical computation is used to calculate these sensitivity functions.

### 4. Simulation example

To illustrate the performances of the proposed OE method for MISO FOS as described in Algorithm 1, a numerical example is presented in this section.

Note that in real implementation, unknown initial conditions can introduce a transient of the simulated model, thus leading to bias on the estimation. Therefore, a history-function-based initialization [13] or alternatively, an infinite-dimensional state-space representation [34] can be used to better initialize fractional systems, so that the parameter estimation gets more consistent.

System identification, and particularly for fractional systems, is usually undertaken only around null initial conditions or around a working point. Therefore, the identified model is valid only for small variations around the considered working point.

To avoid local minima first, a global S-commensurate order is estimated for all subsystems which is used as an initial hit for the second variant where local S-commensurate orders are estimated, which in turn is used as an initial hit for the last variant to estimate all differentiation orders.
The following MISO system is used to generate data:

$$ S' : \left\{ \begin{array}{l} G_1(s) = \frac{1}{1 + 1.5s^{1.2} + 0.5s^{2.8}} \\ G_2(s) = \frac{0.5}{1 + 1.5s^{1.1} + 0.4s^{2.5}} \\ G_3(s) = \frac{1.5}{1 + 1.5s^{0.7} + 0.6s^{1.6}} \end{array} \right. \qquad (4.1) $$

<figure><img src="image_2.png" /><figcaption>FIGURE 2. Input/output data used for the identification of system (4.1)</figcaption></figure>

Three uncorrelated pseudo random binary sequences, $u_i$ for $i = 1, 2, 3$ with amplitudes set between -1 and 1, are applied to the MISO system (4.1). The sampling time is set to $T = 0.2s$. The noise-free output $y(t)$, stemming from the sum of the three outputs, is corrupted by an additive zero mean Gaussian white noise with a noise to signal ratio (NSR) set to -15dB (see Fig. 2 with $K = 3$).

As it is known in system identification, the model structure should be set as the true system in order to get the less modeling error:

$$ G_k = \frac{b_{0,k}}{1 + a_{1,k}s^{\alpha_{1,k}} + a_{2,k}s^{\alpha_{2,k}}}, \quad \text{for } k = 1, 2, 3. \qquad (4.2) $$

All the parameters are to be estimated, namely the differentiation orders and the coefficients of all the subsystems. As any gradient-based methods, the initial guess is always tricky as the OE method converges to the nearest local minimum. Therefore, in order to help the convergence towards the global minimum, it is recommended to proceed with the proposed three stage variants.
The following MISO system is used to generate data:

$$ S' : \left\{ \begin{array}{l} G_1(s) = \frac{1}{1 + 1.5s^{1.2} + 0.5s^{2.8}} \\ G_2(s) = \frac{0.5}{1 + 1.5s^{1.1} + 0.4s^{2.5}} \\ G_3(s) = \frac{1.5}{1 + 1.5s^{0.7} + 0.6s^{1.6}} \end{array} \right. \qquad (4.1) $$

<figure><img src="image_2.png" /><figcaption>FIGURE 2. Input/output data used for the identification of system (4.1)</figcaption></figure>

Three uncorrelated pseudo random binary sequences, $u_i$ for $i = 1, 2, 3$ with amplitudes set between -1 and 1, are applied to the MISO system (4.1). The sampling time is set to $T = 0.2s$. The noise-free output $y(t)$, stemming from the sum of the three outputs, is corrupted by an additive zero mean Gaussian white noise with a noise to signal ratio (NSR) set to -15dB (see Fig. 2 with $K = 3$).

As it is known in system identification, the model structure should be set as the true system in order to get the less modeling error:

$$ G_k = \frac{b_{0,k}}{1 + a_{1,k}s^{\alpha_{1,k}} + a_{2,k}s^{\alpha_{2,k}}}, \quad \text{for } k = 1, 2, 3. \qquad (4.2) $$

All the parameters are to be estimated, namely the differentiation orders and the coefficients of all the subsystems. As any gradient-based methods, the initial guess is always tricky as the OE method converges to the nearest local minimum. Therefore, in order to help the convergence towards the global minimum, it is recommended to proceed with the proposed three stage variants.
<figure><img src="image_3.png" /><figcaption>FIGURE 3. All differentiation order estimation versus number of iterations</figcaption></figure>

The differentiation order evolutions are plotted in Fig. 3: first, a global S-commensurate order $\nu$ (—) is estimated; then local S-commensurate orders are estimated ($\nu_1$ for $G_1$ (—), $\nu_2$ for $G_2$ (—) and $\nu_3$ for $G_3$ (—-*)); finally, the S-commensurability constraint being released, all differentiation orders are estimated.

Fig. 4 illustrates the criterion evolution with respect to the three stages, and confirms that all parameters, as well the differentiation orders as the coefficients, are correctly estimated with the MISO OE algorithm:

* for the first stage, when a global S-commensurate order is estimated for all subsystems, the cost function $J_{dB} = -3.88\text{dB}$, so the $\ell_2$-norm of the modeling error is around 11.12dB (see criterion $J_1$ in Fig. 4);

* for the second stage when local S-commensurate orders are estimated, the cost function $J_{dB} = -7.25\text{dB}$, so the $\ell_2$-norm of the modeling error has decreased to 7.75dB (see criterion $J_2$ in Fig. 4);

* for the last stage, when all differentiation orders are estimated without any constraint, the cost function $J_{dB} = -14.98\text{dB}$, so the $\ell_2$-norm of the modeling error is around zero dB (see criterion $J_3$ in Fig. 4).

A Monte Carlo simulation with 75 runs has been carried out and the results are provided in table 1. The statistical properties of the algorithm are remarkable: the mean values are close to the true parameters and the standard deviations are very low. The proposed MISO OE algorithm for differentiation order and coefficient estimation is therefore validated in this example as the estimated parameters have converged to the true ones in a
<figure><img src="image_3.png" /><figcaption>FIGURE 3. All differentiation order estimation versus number of iterations</figcaption></figure>

The differentiation order evolutions are plotted in Fig. 3: first, a global S-commensurate order $\nu$ (—) is estimated; then local S-commensurate orders are estimated ($\nu_1$ for $G_1$ (—), $\nu_2$ for $G_2$ (—) and $\nu_3$ for $G_3$ (—-*)); finally, the S-commensurability constraint being released, all differentiation orders are estimated.

Fig. 4 illustrates the criterion evolution with respect to the three stages, and confirms that all parameters, as well the differentiation orders as the coefficients, are correctly estimated with the MISO OE algorithm:

* for the first stage, when a global S-commensurate order is estimated for all subsystems, the cost function $J_{dB} = -3.88\text{dB}$, so the $\ell_2$-norm of the modeling error is around 11.12dB (see criterion $J_1$ in Fig. 4);

* for the second stage when local S-commensurate orders are estimated, the cost function $J_{dB} = -7.25\text{dB}$, so the $\ell_2$-norm of the modeling error has decreased to 7.75dB (see criterion $J_2$ in Fig. 4);

* for the last stage, when all differentiation orders are estimated without any constraint, the cost function $J_{dB} = -14.98\text{dB}$, so the $\ell_2$-norm of the modeling error is around zero dB (see criterion $J_3$ in Fig. 4).

A Monte Carlo simulation with 75 runs has been carried out and the results are provided in table 1. The statistical properties of the algorithm are remarkable: the mean values are close to the true parameters and the standard deviations are very low. The proposed MISO OE algorithm for differentiation order and coefficient estimation is therefore validated in this example as the estimated parameters have converged to the true ones in a
<table><caption>TABLE 1. Monte Carlo simulation with 75 runs for the estimation of differentiation orders and coefficients ($\bar{\theta}$ is the mean and $\hat{\sigma}(\theta)$ the standard deviation of the estimates)</caption><thead><tr><th rowspan="2"></th><th rowspan="2">True<br>θ</th><th colspan="2">MISO OE</th></tr><tr><th>θ̃</th><th>σ̂<sub>θ</sub></th></tr></thead><tbody><tr><td>b<sub>0,1</sub></td><td>1</td><td>1.0017</td><td>0.0115</td></tr><tr><td>a<sub>1,1</sub></td><td>1.5</td><td>1.5018</td><td>0.0194</td></tr><tr><td>a<sub>2,1</sub></td><td>0.5</td><td>0.5023</td><td>0.0110</td></tr><tr><td>μ<sub>1,1</sub></td><td>1.2</td><td>1.1979</td><td>0.0117</td></tr><tr><td>μ<sub>2,1</sub></td><td>2.8</td><td>2.7979</td><td>0.0162</td></tr><tr><td>b<sub>0,2</sub></td><td>0.5</td><td>0.4991</td><td>0.0158</td></tr><tr><td>a<sub>1,2</sub></td><td>1.5</td><td>1.4994</td><td>0.0570</td></tr><tr><td>a<sub>2,2</sub></td><td>0.4</td><td>0.4043</td><td>0.0475</td></tr><tr><td>μ<sub>1,2</sub></td><td>1.1</td><td>1.1013</td><td>0.0411</td></tr><tr><td>μ<sub>2,2</sub></td><td>2.5</td><td>2.5087</td><td>0.0979</td></tr><tr><td>b<sub>0,3</sub></td><td>1.5</td><td>1.5032</td><td>0.0374</td></tr><tr><td>a<sub>1,3</sub></td><td>1.5</td><td>1.4844</td><td>0.1022</td></tr><tr><td>a<sub>2,3</sub></td><td>0.6</td><td>0.6228</td><td>0.1276</td></tr><tr><td>μ<sub>1,3</sub></td><td>0.7</td><td>0.6926</td><td>0.0490</td></tr><tr><td>μ<sub>2,3</sub></td><td>1.6</td><td>1.5932</td><td>0.1104</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>FIGURE 4. Cost function versus algorithm iterations</figcaption></figure>

reduced number of iterations (Fig. 4)), and all runs have converged to the true parameters.
<table><caption>TABLE 1. Monte Carlo simulation with 75 runs for the estimation of differentiation orders and coefficients ($\bar{\theta}$ is the mean and $\hat{\sigma}(\theta)$ the standard deviation of the estimates)</caption><thead><tr><th rowspan="2"></th><th rowspan="2">True<br>θ</th><th colspan="2">MISO OE</th></tr><tr><th>θ̃</th><th>σ̂<sub>θ</sub></th></tr></thead><tbody><tr><td>b<sub>0,1</sub></td><td>1</td><td>1.0017</td><td>0.0115</td></tr><tr><td>a<sub>1,1</sub></td><td>1.5</td><td>1.5018</td><td>0.0194</td></tr><tr><td>a<sub>2,1</sub></td><td>0.5</td><td>0.5023</td><td>0.0110</td></tr><tr><td>μ<sub>1,1</sub></td><td>1.2</td><td>1.1979</td><td>0.0117</td></tr><tr><td>μ<sub>2,1</sub></td><td>2.8</td><td>2.7979</td><td>0.0162</td></tr><tr><td>b<sub>0,2</sub></td><td>0.5</td><td>0.4991</td><td>0.0158</td></tr><tr><td>a<sub>1,2</sub></td><td>1.5</td><td>1.4994</td><td>0.0570</td></tr><tr><td>a<sub>2,2</sub></td><td>0.4</td><td>0.4043</td><td>0.0475</td></tr><tr><td>μ<sub>1,2</sub></td><td>1.1</td><td>1.1013</td><td>0.0411</td></tr><tr><td>μ<sub>2,2</sub></td><td>2.5</td><td>2.5087</td><td>0.0979</td></tr><tr><td>b<sub>0,3</sub></td><td>1.5</td><td>1.5032</td><td>0.0374</td></tr><tr><td>a<sub>1,3</sub></td><td>1.5</td><td>1.4844</td><td>0.1022</td></tr><tr><td>a<sub>2,3</sub></td><td>0.6</td><td>0.6228</td><td>0.1276</td></tr><tr><td>μ<sub>1,3</sub></td><td>0.7</td><td>0.6926</td><td>0.0490</td></tr><tr><td>μ<sub>2,3</sub></td><td>1.6</td><td>1.5932</td><td>0.1104</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>FIGURE 4. Cost function versus algorithm iterations</figcaption></figure>

reduced number of iterations (Fig. 4)), and all runs have converged to the true parameters.
## 5. Conclusion

This paper has presented an output-error method for estimating all the parameters of a fractional MISO system with noisy output. Without prior knowledge of the differentiation orders, the proposed MISO OE algorithm estimates the parameters according to the user's will: the coefficients with a global S-commensurate order for the whole system, or the coefficients with local S-commensurate orders (one for each subsystem), or even the coefficients with all differentiation orders without commensurability constraint. In the latter variant, an initialization procedure is proposed, consisting of estimating a global S-commensurate model first, then local S-commensurate orders, and finally all differentiation orders. The proposed algorithms have been validated by Monte Carlo simulation analysis. For future works, it would be interesting to help fix the number of parameters when the MISO model structure is unknown. Also, extensions to cope with colored output noise could be investigated by using hybrid Box-Jenkins MISO models. A further extension of the MISO-oosrivcf algorithm would be tackle parametric identification with time-delay.

## References

[1] T. Abuaisha and J. Kertzscher, Fractional-order modelling and parameter identification of electrical coils. *Fractional Calculus & Applied Analysis* 22, No 1 (2019), 193–216; DOI: 10.1515/fca-2019-0013.

[2] R. Ashurov and S. Umarov, Determination of the order of fractional derivative for subdiffusion equations. *Fractional Calculus & Applied Analysis* 23, No 6 (2020), 1647–1662; DOI: 10.1515/fca-2020-0081.

[3] Z. Belkhatir and T. Laleg-Kirati, Parameters and fractional differentiation orders estimation for linear continuous-time non-commensurate fractional order systems. *Systems & Control Letters* 115 (2018), 26–33; DOI:10.1016/j.sysconle.2018.02.012.

[4] R. Boudarel, J. Delmas and P. Guichet, *Commande optimale des processus*. Tome 2, chapter 11, Dunod (1967), 201-221.

[5] M. Chetoui, M. Thomassin, R. Malti, M. Aoun, S. Najar, M. Abdelkrim and A. Oustaloup, New consistent methods for order and coefficient estimation of continuous-time errors-in-variables fractional models. *Computers & Mathematics with Applications* 66, No 5 (2013), 860–872, DOI: 10.1016/j.camwa.2013.04.028.

[6] J.-F. Duhé, S. Victor, and P. Melchior, Contributions on artificial potential field method for effective obstacle avoidance. *Fractional Calculus & Applied Analysis* 24, No 2 (2021), 421–446; DOI: 10.1515/fca-2021-0019.

[7] P. Eykhoff, *System identification*. John Wiley and Sons (1974).
## 5. Conclusion

This paper has presented an output-error method for estimating all the parameters of a fractional MISO system with noisy output. Without prior knowledge of the differentiation orders, the proposed MISO OE algorithm estimates the parameters according to the user's will: the coefficients with a global S-commensurate order for the whole system, or the coefficients with local S-commensurate orders (one for each subsystem), or even the coefficients with all differentiation orders without commensurability constraint. In the latter variant, an initialization procedure is proposed, consisting of estimating a global S-commensurate model first, then local S-commensurate orders, and finally all differentiation orders. The proposed algorithms have been validated by Monte Carlo simulation analysis. For future works, it would be interesting to help fix the number of parameters when the MISO model structure is unknown. Also, extensions to cope with colored output noise could be investigated by using hybrid Box-Jenkins MISO models. A further extension of the MISO-oosrivcf algorithm would be tackle parametric identification with time-delay.

## References

[1] T. Abuaisha and J. Kertzscher, Fractional-order modelling and parameter identification of electrical coils. *Fractional Calculus & Applied Analysis* **22**, No 1 (2019), 193–216; DOI: 10.1515/fca-2019-0013.

[2] R. Ashurov and S. Umarov, Determination of the order of fractional derivative for subdiffusion equations. *Fractional Calculus & Applied Analysis* **23**, No 6 (2020), 1647–1662; DOI: 10.1515/fca-2020-0081.

[3] Z. Belkhatir and T. Laleg-Kirati, Parameters and fractional differentiation orders estimation for linear continuous-time non-commensurate fractional order systems. *Systems & Control Letters* **115** (2018), 26–33; DOI:10.1016/j.sysconle.2018.02.012.

[4] R. Boudarel, J. Delmas and P. Guichet, *Commande optimale des processus*. Tome 2, chapter 11, Dunod (1967), 201-221.

[5] M. Chetoui, M. Thomassin, R. Malti, M. Aoun, S. Najar, M. Abdelkrim and A. Oustaloup, New consistent methods for order and coefficient estimation of continuous-time errors-in-variables fractional models. *Computers & Mathematics with Applications* **66**, No 5 (2013), 860 – 872, DOI: 10.1016/j.camwa.2013.04.028.

[6] J.-F. Duhé, S. Victor, and P. Melchior, Contributions on artificial potential field method for effective obstacle avoidance. *Fractional Calculus & Applied Analysis* **24**, No 2 (2021), 421–446; DOI: 10.1515/fca-2021-0019.

[7] P. Eykhoff, *System identification*. John Wiley and Sons (1974).
[8] F. Gao, X.J. Lee, F.X. Fei, H.Q. Tong, Y. F. Deng and H. L. Zhao, Identification time-delayed fractional order chaos with functional extrema model via differential evolution. *Expert Systems with Applications* **41**, No 4 (2014), 1601–1608; DOI: 10.1016/j.eswa.2013.08.057.

[9] Z. Gao, X. Lin, and Y. Zheng, System identification with measurement noise compensation based on polynomial modulating function for fractional-order systems with a known time-delay. *ISA transactions* **79** (2018), 62–72; DOI: 10.1016/j.isatra.2018.04.012.

[10] H. Garnier, M. Mensler, and A. Richard, Continuous-time model identification from sampled data. Implementation issues and performance evaluation. *International Journal of Control* **76**, No 13 (2003), 1337–1357; DOI: 10.1080/0020717031000149636.

[11] H. Garnier, L. Wang and P. Young, Direct identification of continuous-time models from sampled data: Issues, basic solutions and relevance. In: *Identification of continuous-time models from sampled data*. Springer (2008) pp. 1–29.

[12] G. Goodwin and R. Payne, *Dynamic system identification. Experiment design and data analysis*. Academic Press, (1977).

[13] T. Hartley, C. Lorenzo, J.-C. Trigeassou, and N. Maamri, Equivalence of history-function based and infinite-dimensional-state initializations for fractional-order operators. *ASME Journal of computational and nonlinear dynamics* **8** (2013) No. 4, 041014–041014–7.

[14] T. Higo, T. Kawaguchi, M. Itamoto and S. Adachi, Recursive identification of fractional time-delay systems using discrete-time model. *IFAC-PapersOnLine* **51 15** (2018), 939–944.

[15] M. Ichise, Y. Nagayanagi and T. Kojima, An analog simulation of non-integer order transfer functions for analysis of electrode processes. *Journal of Electroanalytical Chemistry and Interfacial Electrochemistry* **33**, No 2 (1971), 253–265; DOI: 10.1016/S0022-0728(71)80115-8.

[16] M.N. Kapetina, M.R. Rapaić, A. Pisano and Z.D. Jeličić, Adaptive parameter estimation in LTI systems. *IEEE Transactions on Automatic Control* **64**, No 10 (2019), 4188–4195; DOI: 10.1109/TAC.2019.2893973.

[17] G. Liu, L. Wang, W. Luo, J. Liu and Z. Lu, Parameter identification of fractional order system using enhanced response sensitivity approach. *Communications in Nonlinear Science and Numerical Simulation* **67** (2019), 492–505; DOI: 10.1016/j.cnsns.2018.07.026.

[18] L. Ljung, *System identification – Theory for the user*, 2 ed. Prentice-Hall, Upper Saddle River, N.J., USA (1999).
[8] F. Gao, X.J. Lee, F.X. Fei, H.Q. Tong, Y. F. Deng and H. L. Zhao, Identification time-delayed fractional order chaos with functional extrema model via differential evolution. *Expert Systems with Applications* **41**, No 4 (2014), 1601–1608; DOI: 10.1016/j.eswa.2013.08.057.

[9] Z. Gao, X. Lin, and Y. Zheng, System identification with measurement noise compensation based on polynomial modulating function for fractional-order systems with a known time-delay. *ISA transactions* **79** (2018), 62–72; DOI: 10.1016/j.isatra.2018.04.012.

[10] H. Garnier, M. Mensler, and A. Richard, Continuous-time model identification from sampled data. Implementation issues and performance evaluation. *International Journal of Control* **76**, No 13 (2003), 1337–1357; DOI: 10.1080/0020717031000149636.

[11] H. Garnier, L. Wang and P. Young, Direct identification of continuous-time models from sampled data: Issues, basic solutions and relevance. In: *Identification of continuous-time models from sampled data*. Springer (2008) pp. 1–29.

[12] G. Goodwin and R. Payne, *Dynamic system identification. Experiment design and data analysis*. Academic Press, (1977).

[13] T. Hartley, C. Lorenzo, J.-C. Trigeassou, and N. Maamri, Equivalence of history-function based and infinite-dimensional-state initializations for fractional-order operators. *ASME Journal of computational and nonlinear dynamics* **8** (2013) No. 4, 041014–041014–7.

[14] T. Higo, T. Kawaguchi, M. Itamoto and S. Adachi, Recursive identification of fractional time-delay systems using discrete-time model. *IFAC-PapersOnLine* **51 15** (2018), 939–944.

[15] M. Ichise, Y. Nagayanagi and T. Kojima, An analog simulation of non-integer order transfer functions for analysis of electrode processes. *Journal of Electroanalytical Chemistry and Interfacial Electrochemistry* **33**, No 2 (1971), 253–265; DOI: 10.1016/S0022-0728(71)80115-8.

[16] M.N. Kapetina, M.R. Rapaić, A. Pisano and Z.D. Jeličić, Adaptive parameter estimation in LTI systems. *IEEE Transactions on Automatic Control* **64**, No 10 (2019), 4188–4195; DOI: 10.1109/TAC.2019.2893973.

[17] G. Liu, L. Wang, W. Luo, J. Liu and Z. Lu, Parameter identification of fractional order system using enhanced response sensitivity approach. *Communications in Nonlinear Science and Numerical Simulation* **67** (2019), 492–505; DOI: 10.1016/j.cnsns.2018.07.026.

[18] L. Ljung, *System identification – Theory for the user*, 2 ed. Prentice-Hall, Upper Saddle River, N.J., USA (1999).
[19] R.L. Magin Fractional calculus models of complex dynamics in bi-
ological tissues. *Comput. Math. Appl.* **59** (2010), 1586–1593; DOI:
10.1016/j.camwa.2009.08.039.

[20] R. Mansouri, M. Bettayeb and S. Djennoune, Multivariable fractional
system approximation with initial conditions using integral state space
representation. *Computers & Mathematics with Applications* **59**, No 5
(2010), 1842 – 1851; DOI: 10.1016/j.camwa.2009.08.024.

[21] D. Matignon, Stability properties for generalized fractional differential
systems. *ESAIM proceedings - Systèmes Différentiels Fractionnaires -
Modèles, Méthodes et Applications* **5** (1998).

[22] A. Narang, S. Shah and T. Chen, Continuous-time model identification of fractional-order models with time delays. *IET Control Theory & Applications* **5**, No 7 (2011), 900–912; DOI: 10.1049/iet-cta.2010.0718.

[23] I. Podlubny, *Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications*. Academic Press, San Diego (1999).

[24] T. Poinot and J.-C. Trigeassou, Identification of fractional systems using an output-error technique. *Nonlinear Dynamics* **38**, No 1 (2004), 133–154; DOI: 10.1007/s11071-004-3751-y.

[25] G. Rao and H. Unbehauen, Identification of continuous-time systems.
*Electrical Engineering* **153** (2006), 185–220; DOI:.

[26] M.R. Rapaić and R. Malti, On stability regions of fractional systems
in the space of perturbed orders. *IET Control Theory & Applications*
**13**, No 16 (2019); DOI: 10.1049/iet-cta.2018.6350.

[27] M.R. Rapaić and A. Pisano, Variable-order fractional operators for adaptive order and parameter estimation. *IEEE Transactions on Automatic Control* **59**, No 3 (March 2014), 798–803; DOI: 10.1109/TAC.2013.2278136.

[28] J.-B. Receveur, S. Victor and P. Melchior, New interpretation of frac-
tional potential fields for robust path planning. *Fractional Calculus &
Applied Analysis* **22**, No 1 (2019), 113–127; DOI: 10.1515/fca-2019-
0007.

[29] J. Richalet, *Pratique de l'identification*. Hermès (1991)

[30] J. Richalet, A. Rault and R. Pouliquen, *Identification des processus par la méthode du modèle*. Gordon and Breach (1971)

[31] T. Söderström and P. Stoica, *System Identification. Series in Systems and Control Engineering*. Prentice-Hall (1989).

[32] Y. Tang, H. Liu, W. Wang, Q. Lian and X. Guan, Parameter identi-
fication of fractional order systems using block pulse functions. *Signal
Processing* **107**, (2015), 272 – 281; DOI: 10.1016/j.sigpro.2014.04.011.
[19] R.L. Magin Fractional calculus models of complex dynamics in bi-ological tissues. *Comput. Math. Appl.* **59** (2010), 1586–1593; DOI: 10.1016/j.camwa.2009.08.039.

[20] R. Mansouri, M. Bettayeb and S. Djennoune, Multivariable fractional system approximation with initial conditions using integral state space representation. *Computers & Mathematics with Applications* **59**, No 5 (2010), 1842–1851; DOI: 10.1016/j.camwa.2009.08.024.

[21] D. Matignon, Stability properties for generalized fractional differential systems. *ESAIM proceedings - Systèmes Différentiels Fractionnaires - Modèles, Méthodes et Applications* **5** (1998).

[22] A. Narang, S. Shah and T. Chen, Continuous-time model identification of fractional-order models with time delays. *IET Control Theory & Applications* **5**, No 7 (2011), 900–912; DOI: 10.1049/iet-cta.2010.0718.

[23] I. Podlubny, *Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications*. Academic Press, San Diego (1999).

[24] T. Poinot and J.-C. Trigeassou, Identification of fractional systems using an output-error technique. *Nonlinear Dynamics* **38**, No 1 (2004), 133–154; DOI: 10.1007/s11071-004-3751-y.

[25] G. Rao and H. Unbehauen, Identification of continuous-time systems. *Electrical Engineering* **153** (2006), 185–220; DOI:.

[26] M.R. Rapaić and R. Malti, On stability regions of fractional systems in the space of perturbed orders. *IET Control Theory & Applications* **13**, No 16 (2019); DOI: 10.1049/iet-cta.2018.6350.

[27] M.R. Rapaić and A. Pisano, Variable-order fractional operators for adaptive order and parameter estimation. *IEEE Transactions on Automatic Control* **59**, No 3 (March 2014), 798–803; DOI: 10.1109/TAC.2013.2278136.

[28] J.-B. Receveur, S. Victor and P. Melchior, New interpretation of frac-tional potential fields for robust path planning. *Fractional Calculus & Applied Analysis* **22**, No 1 (2019), 113–127; DOI: 10.1515/fca-2019-0007.

[29] J. Richalet, *Pratique de l'identification*. Hermès (1991)

[30] J. Richalet, A. Rault and R. Pouliquen, *Identification des processus par la méthode du modèle*. Gordon and Breach (1971)

[31] T. Söderström and P. Stoica, *System Identification. Series in Systems and Control Engineering*. Prentice-Hall (1989).

[32] Y. Tang, H. Liu, W. Wang, Q. Lian and X. Guan, Parameter identi-fication of fractional order systems using block pulse functions. *Signal Processing* **107**, (2015), 272 – 281; DOI: 10.1016/j.sigpro.2014.04.011.
[33] M. Thomassin and R. Malti, Subspace method for continuous-time fractional system identification. *IFAC Proceedings Volumes - 15th IFAC Symposium on System Identification (SYSID'2009)* **42**, No 10 (2009), 880-885; DOI: 10.3182/20090706-3-FR-2004.00146.

[34] J. C. Trigeassou and N. Maamri, Initial conditions and initialization of linear fractional differential equations. *Signal Processing* **91**, No. 3 (2011), 427-436.

[35] S. Victor, R. Malti, H. Garnier and A. Oustaloup, Parameter and differentiation order estimation in fractional models. *Automatica* **49**, No 4 (2013), 926-935; DOI: 10.1016/j.automatica.2013.01.026.

[36] S. Victor and P. Melchior, Improvements on flat output characterization for fractional systems. *Fractional Calculus & applied analysis* **18**, No 1 (2015), 238-260; DOI: 10.1515/fca-2015-0016.

[37] S. Victor, P. Melchior, J. Lévine and A. Oustaloup, Flatness for linear fractional systems with application to a thermal system. *Automatica* **57** (2015), 213-221; DOI: 10.1016/j.automatica.2015.04.021.

[38] S. Victor, P. Melchior, M. Pellet and A. Oustaloup, Lung thermal transfer system identification with fractional models. *IEEE Transactions on Control Systems Technology* **28**, No 1 (2020), 172-182; DOI: 10.1109/TCST.2018.2877606.

[39] B. Wahlberg, Identification of linear systems: A practical guideline to accurate modeling: J. Schoukens and R. Pintelon. *Autom. 29* (1993), 1621-1622; DOI: 10.1016/0005-1098(93)90036-S.

[40] Y. Wang, Y. Chen and X. Liao, State-of-art survey of fractional order modeling and estimation methods for lithium-ion batteries. *Fractional Calculus & Applied Analysis* **22**, No 6 (2019), 1449-1479; DOI: 10.1515/fca-2019-0076.

[41] Y. Wei, D. Liu, D. Boutat and Y. Chen, An improved pseudo-state estimator for a class of commensurate fractional order linear systems based on fractional order modulating functions. *Systems & Control Letters* **118** (2018), 29-34; DOI: 10.1016/j.sysconle.2018.05.011.

[42] P. C. Young, *Recursive estimation and Time-Series Analysis*. Springer-Verlag, Berlin (1984).

[43] P.C. Young and H. Garnier, Identification and estimation of continuous-time data-based mechanistic (dbm) models for environmental systems. *Environmental Modelling and Software* **21**, No 8 (2006), 1055-1072; DOI: 10.1016/j.envsoft.2005.05.007.

<footnote n="1">¹ *Univ. Bordeaux, IMS - UMR 5218 CNRS, Bordeaux INP/ENSEIRB-MATMECA,*
*351 Cours de la Libération, 33405 Talence CEDEX, France*</footnote>
[33] M. Thomassin and R. Malti, Subspace method for continuous-time fractional system identification. *IFAC Proceedings Volumes - 15th IFAC Symposium on System Identification (SYSID'2009)* **42**, No 10 (2009), 880-885; DOI: 10.3182/20090706-3-FR-2004.00146.

[34] J. C. Trigeassou and N. Maamri, Initial conditions and initialization of linear fractional differential equations. *Signal Processing* **91**, No. 3 (2011), 427-436.

[35] S. Victor, R. Malti, H. Garnier and A. Oustaloup, Parameter and differentiation order estimation in fractional models. *Automatica* **49**, No 4 (2013), 926-935; DOI: 10.1016/j.automatica.2013.01.026.

[36] S. Victor and P. Melchior, Improvements on flat output characterization for fractional systems. *Fractional Calculus & applied analysis* **18**, No 1 (2015), 238-260; DOI: 10.1515/fca-2015-0016.

[37] S. Victor, P. Melchior, J. Lévine and A. Oustaloup, Flatness for linear fractional systems with application to a thermal system. *Automatica* **57** (2015), 213-221; DOI: 10.1016/j.automatica.2015.04.021.

[38] S. Victor, P. Melchior, M. Pellet and A. Oustaloup, Lung thermal transfer system identification with fractional models. *IEEE Transactions on Control Systems Technology* **28**, No 1 (2020), 172-182; DOI: 10.1109/TCST.2018.2877606.

[39] B. Wahlberg, Identification of linear systems: A practical guideline to accurate modeling: J. Schoukens and R. Pintelon. *Autom. 29* (1993), 1621-1622; DOI: 10.1016/0005-1098(93)90036-S.

[40] Y. Wang, Y. Chen and X. Liao, State-of-art survey of fractional order modeling and estimation methods for lithium-ion batteries. *Fractional Calculus & Applied Analysis* **22**, No 6 (2019), 1449-1479; DOI: 10.1515/fca-2019-0076.

[41] Y. Wei, D. Liu, D. Boutat and Y. Chen, An improved pseudo-state estimator for a class of commensurate fractional order linear systems based on fractional order modulating functions. *Systems & Control Letters* **118** (2018), 29-34; DOI: 10.1016/j.sysconle.2018.05.011.

[42] P. C. Young, *Recursive estimation and Time-Series Analysis*. Springer-Verlag, Berlin (1984).

[43] P.C. Young and H. Garnier, Identification and estimation of continuous-time data-based mechanistic (dbm) models for environmental systems. *Environmental Modelling and Software* **21**, No 8 (2006), 1055-1072; DOI: 10.1016/j.envsoft.2005.05.007.

<footnote n="1">¹ *Univ. Bordeaux, IMS - UMR 5218 CNRS, Bordeaux INP/ENSEIRB-MATMECA,*
*351 Cours de la Libération, 33405 Talence CEDEX, France*</footnote>
(e-mail: abir.mayoufi@ims-bordeaux.fr, stephane.victor@ims-bordeaux.fr (Corr. author), rachid.malti@ims-bordeaux.fr).

² University of Gabes, ENIG, MACS LR16ES22

Avenue Omar Ibn El Khattab, Zrig, 6029 Gabes, Tunisia

(e-mail: abir.mayoufi@macs.tn, chetoui.manel@macs.tn, mohamed.aoun@macs.tn).
(e-mail: abir.mayoufi@ims-bordeaux.fr, stephane.victor@ims-bordeaux.fr (Corr. author), rachid.malti@ims-bordeaux.fr).

² University of Gabes, ENIG, MACS LR16ES22

Avenue Omar Ibn El Khattab, Zrig, 6029 Gabes, Tunisia

(e-mail: abir.mayoufi@macs.tn, chetoui.manel@macs.tn, mohamed.aoun@macs.tn).