# Contributions on artificial potential field method for effective obstacle avoidance

Jean-François Duhé, Stephane Victor, Pierre Melchior

► To cite this version:

Jean-François Duhé, Stephane Victor, Pierre Melchior. Contributions on artificial potential field method for effective obstacle avoidance. Fractional Calculus and Applied Analysis, 2021, 24 (2), pp.421-446. 10.1515/fca-2021-0019. hal-03548486

HAL Id: hal-03548486

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

Submitted on 9 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
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 num-
ber 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 al-
gorithm. 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 re-
leased, 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 es-
timated 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 mode-
ling diffusive phenomena since the 1950s: thermal engineering [23], elec-
trochemistry [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 frac-
tional (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{l}{h} \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 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:

$$p f(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}{h} \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*.

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), \quad (2.4)$$

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

<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>
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Fractional MISO system structure</figcaption></figure>

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].
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^{3.6 \times 0.4} \cdot 7 \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 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:

<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>
$$ \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 $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.

### 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.

<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}^{iter}$ as

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

with

$$ \left\{ \begin{array}{l} \frac{\partial J}{\partial \hat{\theta}} = \sum_{h=1}^{H} \frac{\partial \varepsilon(t_h)^T}{\partial \hat{\theta}} \varepsilon(t_h) : \text{gradient} \\ \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} \\ \xi : \text{Marquardt parameter} \end{array} \right. \quad (3.5) $$

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

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

**Step 3: 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^{iter}(t) \quad (3.7) $$

with $y_k^{iter}(t)$ being calculated at the last iteration.
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) $$

for $k = 1, \dots, 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, \dots, M_k \qquad (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. \qquad (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}, \qquad (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.

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) $$

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 2. Input/output data used for the identification of system (4.1)</figcaption></figure>

<figure><img src="image_4.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.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Cost function versus algorithm iterations</figcaption></figure>

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.12\text{dB}$ (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.75\text{dB}$ (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 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>True</th><th colspan="2">MISO OE</th></tr><tr><th>θ</th><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>

## 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. *Fract. Calc. Appl. Anal.* **22**, No 1 (2019), 193–216; DOI: 10.1515/fca-2019-0013; https://www.degruyter.com/journal/key/fca/21/1/html.

[2] R. Ashurov and S. Umarov, Determination of the order of fractional derivative for subdiffusion equations. *Fract. Calc. Appl. Anal.* **23**, No 6 (2020), 1647–1662; DOI: 10.1515/fca-2020-0081; https://www.degruyter.com/journal/key/fca/23/6/html.

[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, Ch. 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 & Math. with Appl.* **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. *Fract. Calc. Appl. Anal.* **24**, No 2 (2021), 421–446; DOI: 10.1515/fca-2021-0019; https://www.degruyter.com/journal/key/fca/24/2/html.

[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 Appl.* **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 Trans.* **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. *Intern. J. 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), 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 J. of Comput. and Nonlin. 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. *J. of Electroanal. Chemistry and Interfacial Electrochem.* **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 Trans. 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. *Commun. in Nonlin. Sci. and Numer. Simul.* **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 biological 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 & Math. with Appl.* **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 Proc. - 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 & Appl.* **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. *Nonlin. 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. *Elect. Engineering* **153** (2006), 185-220.

[26] M.R. Rapaić and R. Malti, On stability regions of fractional systems in the space of perturbed orders. *IET Control Theory & Appl.* **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 Trans. on Automatic Control* **59**, No 3 (2014), 798-803; DOI: 10.1109/TAC.2013.2278136.

[28] J.-B. Receveur, S. Victor and P. Melchior, New interpretation of fractional potential fields for robust path planning. *Frac. Calc. Appl. Anal.* **22**, No 1 (2019), 113-127; DOI: 10.1515/fca-2019-0007; https://www.degruyter.com/journal/key/fca/22/1/html.

[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 identification 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 Proc. Volumes - 15th IFAC Symp. 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. *Fract. Calc. Appl. Anal.* **18**, No 1 (2015), 238-260; DOI: 10.1515/fca-2015-0016; https://www.degruyter.com/journal/key/fca/18/1/html.

[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 Trans. 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. *Fract. Calc. Appl. Anal.* **22**, No 6 (2019), 1449–1479; DOI: 10.1515/fca-2019-0076; https://www.degruyter.com/journal/key/fca/22/6/html.

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

Received: June 5, 2021

rachid.malti@ims-bordeaux.fr

<footnote n="1">¹ Univ. Bordeaux, IMS - UMR 5218 CNRS
Bordeaux INP/ENSEIRB-MATMECA,
351 Cours de la Libération, 33405 Talence Cedex, FRANCE
e-mails: abir.mayoufi@ims-bordeaux.fr
stephane.victor@ims-bordeaux.fr (Corresponding author)</footnote>

<footnote n="2">² University of Gabes, ENIG, MACS LR16ES22
Avenue Omar Ibn El Khattab, Zrig, 6029 Gabes, TUNISIA
e-mails: abir.mayoufi@macs.tn,
chetoui.manel@macs.tn, mohamed.aoun@macs.tn</footnote>