# A comparison between two methods for MISO fractional models estimation

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

► To cite this version:

Abir Mayoufi, Manel Chetoui, Stéphane Victor, Mohamed Aoun, Rachid Malti. A comparison between two methods for MISO fractional models estimation. 2020 17th International Multi-Conference on Systems, Signals & Devices (SSD), Jul 2020, Monastir, Tunisia. pp.446-451, 10.1109/SSD49366.2020.9364187. hal-03173085

HAL Id: hal-03173085

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

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
# A comparison between two methods for MISO fractional models estimation

Abir MAYOUFI

Research Laboratory Modeling, Analysis and Control of Systems (LR16ES22)
National Engineering School of Gabes
University of Gabes, Tunisia
abir.maayoufi@gmail.com

Manel CHETOUI

Research Laboratory Modeling, Analysis and Control of Systems (LR16ES22)
National Engineering School of Gabes
University of Gabes, Tunisia
manel.chetoui@gmail.com

Stéphane VICTOR

Laboratoire IMS UMR 5218 CNRS
Equipe CRONE
Université de Bordeaux
F-33405 Talence, France
stephane.victor@u-bordeaux.fr

Mohamed AOUN

Research Laboratory Modeling, Analysis and Control of Systems (LR16ES22)
National Engineering School of Gabes
University of Gabes, Tunisia
mohamed.aoun@gmail.com

Rachid MALTI

Laboratoire IMS UMR 5218 CNRS
Equipe CRONE
Université de Bordeaux
F-33405 Talence, France
rachid.malti@u-bordeaux.fr

**Abstract**—This paper proposes two new methods for multiple input-single output system identification with fractional models: the instrumental variable based method and the output-error based method. The fractional orders are supposed known and the linear coefficients are estimated. A comparative study between the developed methods is illustrated via a numerical example. Monte Carlo simulations are used to demonstrate the efficiency of the two methods.

***Index Terms***—system identification, continuous-time, output-error, instrumental variable, multiple input single output model, fractional derivative, state variable filter.

## I. INTRODUCTION

The objective of modelling is to approximate the system from the data measured at the output and/or input. The used models then describe a relationship between the input and the output of the system in order to predict the output of the system at the next moment. This relationship is defined through parameters that must be determined for each system. This step of calculating the parameters is the system identification.

Among the methods of system identification, there is the **Instrumental Variable (IV)** method which was introduced for the first time by [1] in Discrete-Time (DT). In [2] this method was extended to the Continuous-Time case, all this work being limited for single-input-single-output (**SISO**) systems.

First extension of the **Refined Instrumental Variables** for **Continuous-time** system method to the multiple-input-single-output system was made by [5].

Another method of system identification can be mentioned based on the **output-error method**. In [6] the identification of output error models from frequency domain data generally results in a non-convex optimization

problem was proposed, in [17] an output error method is proposed for the identification of continuous-time systems with time delay, in this work the **output-error method** was introduced in both single-input-single-output system and multiple-input-multiple-output system **MIMO**.

The use of the fractional differentiation operator is more and more widespread in fields such as mechanics, automation, thermal or electrochemistry, the late nineties was the beginning of fractional system identification and continues until now. Among fractional system identification works, there are [7] [8] [15] [14]. This was not limited by those research and it was a profound by other publication [11] [12] [13] [16].

In [9] the **Instrumental Variables** for **Continuous-time** method was extended to the single-input-single-output fractional system since it provided consistent estimates in presence of white additive noise.

In this paper an extension of the both methods to handel **MISO** fractional model identification is suggested. With a comparison between the two methods to demonstrate the performances of each method and clarify the most efficient.

The paper is organized as follows. After introducing the paper in section 1, fractional calculus was introduced in section 2, the problem formulation is presented in section 3. In section 4 and 5 the algorithm of each approach is detailed. Section 6 validates the performances of the developed methods through a numerical example. Finally a conclusion is drawn in section 7.

## II. MATHEMATICAL BACKGROUND

A fractional mathematical **SISO** model can be described by a differential equation containing operators of
fractional order:

$$y(t) + a_1 p^{\alpha_1} y(t) + \dots + a_N p^{\alpha_N} y(t) = b_0 p^{\beta_0} u(t) + b_1 p^{\beta_1} u(t) + \dots + b_M p^{\beta_M} u(t) \quad (1)$$

where $u(t)$, $y(t)$ respectively are the input and the noise-free output, $(a_i(i = 1, \dots, N), b_j(j = 0, \dots, M)) \in \mathbb{R}^2$ and the differentiation orders are ordered positive real numbers such as:

$$\begin{cases} 0 < \alpha_1 < \alpha_2 < \dots < \alpha_N \\ 0 < \beta_0 < \beta_1 < \dots < \beta_M \end{cases}$$

$p$ refers to the differential operator ($p = \frac{d}{dt}$)).

The concept of differentiation to an arbitrary order $\nu$ of a function $f(t)$, in the sense of Grünwald-Letnikov, is defined by:

$$p^{\nu} f(t) = \frac{1}{T^{\nu}} \sum_{h=0}^{H} (-1)^h \begin{pmatrix} h \\ \nu \end{pmatrix} f(t - Th) \quad (2)$$

with $H = \lfloor \frac{t}{T} \rfloor$ ($\lfloor \cdot \rfloor$ is the floor operator), $T$ is the sampling time and $\begin{pmatrix} \nu \\ k \end{pmatrix}$ is the Newton's binomial generalized to the fractional orders:

$$\begin{pmatrix} \nu \\ k \end{pmatrix} = \frac{\Gamma(\nu + 1)}{\Gamma(k + 1)\Gamma(\nu - k + 1)} = \frac{\nu(\nu - 1)\dots(\nu - k + 1)}{k!}$$

The Laplace transform of a $\nu$ order derivative, of a causal function $f$ (null $\forall t \le 0$), is given by:

$$\mathcal{L}\{p^{\nu}x(t)\} = s^{\nu}X(s) \quad (3)$$

$s$ is the Laplace variable.

This leads to the following non-integer transfer function:

$$G(s) = \frac{B(s)}{A(s)} = \frac{\sum_{j=0}^{M} b_j s^{\beta_j}}{1 + \sum_{i=1}^{N} a_i s^{\alpha_i}} \quad (4)$$

For a commensurate order system, the transfer function can be rewritten in the following form:

$$G(s) = \frac{\sum_{j=0}^{m} \tilde{b}_j s^{j\nu}}{1 + \sum_{i=1}^{n} \tilde{a}_i s^{i\nu}}; \nu \in \mathbb{R}^* \quad (5)$$

with $j = 1, \dots, m$; $m = \frac{\beta_M}{\nu} \in \mathbb{N}$ and $i = 1, \dots, n$; $n = \frac{\alpha_N}{\nu} \in \mathbb{N}$.

$$\left\{ \begin{array}{ll} \tilde{b}_j = b_j & \text{if } \exists j \in \{0, 1, \dots, M\} \text{ such that } j\nu = \beta_j \\ \tilde{b}_j = 0 & \text{otherwise} \\ \tilde{a}_i = a_i & \text{if } \exists i \in \{1, \dots, N\} \text{ such that } i\nu = \alpha_i \\ \tilde{a}_i = 0 & \text{otherwise} \end{array} \right.$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Fractional MISO model.</figcaption></figure>

### III. PROBLEM FORMULATION

Consider a MISO fractional system (Fig 1), represented by:

$$\left\{ \begin{aligned} y_{u_k}(t) &= G_k(p)u_k(t) \\ y(t) &= \sum_{k=1}^{n_u} y_{u_k}(t) \\ y^*(t_h) &= y(t_h) + e(t_h) \end{aligned} \right. \quad (6)$$

where $n_u$ is the number of inputs and $u(t) = [u_1(t), \dots, u_{n_u}(t)]$. The input-output data set $(u(t), y(t))$ are collected at regular samples. The input signals $u_1(t), u_2(t), \dots, u_{n_u}(t)$ are supposed to be uncorrelated. The output measurements $y(t)$ are affected by an additive white noise $e(t_h)$ and $G_k(p)$ is the fractional transfer function that relates the input signal $u_k(t)$ to its noise-free output $y_{u_k}(t)$, and defined as follows:

$$G_k(p) = \frac{B_k(p, \theta_k)}{A_k(p, \theta_k)} = \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}}} \quad (7)$$

The objective is to estimate the model parameter represented by equation (6) using $H$ samples of inputs and output data $\{u_1(t_h), \dots, u_{n_u}(t_h), y^*(t_h)\}_{h=1}^H$. The fractional commensurate order is supposed known a priori and only the linear coefficients are estimated. The parameter vector is defined as:

$$\left\{ \begin{array}{l} \theta_k = [b_{0,k}, b_{1,k}, \dots, b_{M_k,k}, a_{1,k}, \dots, a_{N_k,k}]^T, \quad k = 1, \dots, n_u \\ \theta = [\theta_1; \dots; \theta_{n_u}] \end{array} \right. \quad (8)$$

It is well known that the methods based on least squares give biased parameters in the presence of important noise affecting the output measurements. To solve this problem, the use of methods based on the output-error and the instrumental variable are proposed in this work.

### IV. OUTPUT ERROR METHOD FOR CONTINUOUS-TIME FRACTIONAL MODELS

A methodology for estimating the linear coefficients based on an output-error method is proposed. This algorithm being iterative, the parameter vector $\theta$ expressed in equation (8) is denoted by $\theta^{iter}$.
The summary of the oe is described in the following steps:

### Step1: Initialization

Use an initial fractional model arbitrary to generate an initial estimated parameters vector $θ^{iter}$ (iter = 1).

### Step2: Iterative Levenberg-Marquardt algorithm

a) Minimize the least-squares curve fitting problem using the initial parameters vector $θ^{iter}$.

b) for iter = 2 to convergence

i) Run the optimisation algorithm Levenberg-Marquardt which interpolates between the Gauss-Newton algorithm and the method of gradient descent, by minimizing the sum of square deviations to converges to a local minimum.

ii) Evaluate the output error at each iteration.

$$\varepsilon(t) = y^*(t) - \sum_{k=1}^{n_u} y_{u_k}^{iter}(t) \quad (9)$$

### Step3: Compute the error covariance matrix using the approximated Hessian.

for $k = 1, ..., n_u$

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

with $\hat{\sigma}^2$ is the empirical estimation of the noise variance, given by:

$$\hat{\sigma}^2 = y^*(t) - \sum_{k=1}^{n_u} y_{u_k}^{iter}(t) \quad (11)$$

where $y_{u_k}^{iter}(t)$ is calculated at the last iteration.

$\mathcal{H}_k$ is the approximated Hessian given by:

$$\mathcal{H}_k = \frac{\partial \varepsilon}{\partial \theta_k} \quad \frac{\partial \varepsilon}{\partial \theta_k} \quad (12)$$

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

$$\frac{\partial \varepsilon}{\partial \theta_k} = -\frac{\partial \hat{y}_{u_k}}{\partial \theta_k} = - \left[ \frac{\partial \hat{y}_{u_k}}{\partial b_{0,k}}, \dots, \frac{\partial \hat{y}_{u_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]^T \quad (13)$$

$\frac{\partial \hat{y}_{u_k}}{\partial b_{j,k}}$ and $\frac{\partial \hat{y}_{u_k}}{\partial a_{i,k}}$ denote the output sensitivity functions, given by:

$$\frac{\partial \hat{y}_{u_k}}{\partial b_{j,k}} = \frac{\hat{p}^{\beta_j}}{1 + \sum_{i=1}^{N_k} \hat{a}_{i,k} \hat{p}^{\hat{\alpha}_i}} u_k(t), \quad \forall j = 0, \dots, M_k \quad (14)$$

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

## V. INSTRUMENTAL VARIABLE STATE-VARIABLE FOR CONTINUOUS-TIME FRACTIONAL MODELS

### A. ivsvf method for SISO models

The ivsvf method is developed in [4] for SISO (assuming that $n_u = 1$ in Fig.1) system identification by

fractional models. It is based on the use of an auxiliary model.

The error function $\varepsilon$, that yields to an optimal instrumental variable estimates as mentioned in [3], is given by the output error equation:

$$\varepsilon(t, \theta) = y^*(t) - \frac{B(p, \theta)}{A(p, \theta)} u(t) \quad (16)$$

Equation (16) can be rewritten in the following form,

$$\varepsilon(t, \theta) = A(p, \theta) \left( \frac{1}{A(p, \theta)} y^*(t) \right) - B(p, \theta) \left( \frac{1}{A(p, \theta)} u(t) \right) \quad (17)$$

The proposed idea in [10] is to use a filtered signals of the input and output instead of the measured ones by applying the following filter:

$$F(p) = \frac{1}{A(p, \theta)} \quad (18)$$

But the coefficients of this filter are unknown. To solve this problem, the use of the state variable filter (svf) has been proposed in [10].

The svf transfer function is defined as:

$$F_{svf}(p) = \frac{1}{\left(1 + \frac{p}{\omega_c}\right)^{N_c}} \quad (19)$$

$\omega_c$ denotes the cut-off frequency and $N_c$ is the filter integer order defined as: $N_c = \lfloor \alpha_N \rfloor + 1$

Applying the svf filter on both input and output signals leads to the filtered signals:

$$\left\{ \begin{array}{l} y_f^*(t) = F_{svf}(p)y^*(t) \\ u_f(t) = F_{svf}(p)u(t) \end{array} \right. \quad (20)$$

In this case, equation (17) takes the following form:

$$\varepsilon(t, \theta) = A(p, \theta)y_f^*(t) + B(p, \theta)u_f(t) \quad (21)$$

The algorithm of the ivsvf method is extended in the next section for MISO system identification with fractional models. The different steps of the algorithm are given for $k = 1$.

### B. ivsvf method for MISO models

The main idea is to decompose the MISO model, described in equation (6), into $n_u$ SISO models. In this case, the output error function takes the following form:

$$\varepsilon_k(t, \theta) = x_{k,f}(t, \theta) - y_{u_k,f}(t, \theta_k), \quad k = 1, \dots, n_u \quad (22)$$

where, $y_{u_k}$ is the noise-free part of $x_k$:

$$x_k(t, \theta) = y^*(t) - \sum_{\substack{n=1 \\ n \neq k}}^{n_u} y_{u_n}(t, \theta_n) \quad (23)$$

noting that:

$$x_{k,f}(t, \theta) = F_{svf}(p)x_k(t, \theta) \quad (24)$$
At each iteration a “for” loop is launched to estimate $\theta_k$ ($k = 1, ..., n_u$), assuming that $\theta_n$ (with $n \neq k$) is known.

In this case, equation (22) can be expressed in the regression form:

$$\varepsilon_k(t, \theta) = x_{k,f}(t, \theta_k) - \varphi_{k,f}(t)\theta_k, \quad k = 1, ..., n_u \quad (25)$$

where,

$$\varphi_{k,f}(t, \theta_k) = \begin{bmatrix} p^{\beta_0} u_{k,f}(t) \dots p^{\beta_M} u_{k,f}(t) \\ -p^{\alpha_1} x_{k,f}(t, \theta_k) \dots -p^{\alpha_N} x_{k,f}(t, \theta_k) \end{bmatrix}^T \quad (26)$$

The MISO version of ivsf method is summarized by the following steps.

### Step1: Initialization

Initialize the parameters vector $\theta^{iter}$ ($iter = 0$) to calculate the auxiliary models outputs $y_{u_k}$.

$$y_{u_k}(t, \theta_k^{iter}) = \frac{B_k(p, \theta_k^{iter})}{A_k(p, \theta_k^{iter})} u_k(t), \quad k = 1, ..., n_u \quad (27)$$

### Step2:

$iter = 1, ..., N_{iter}$

A) for $k = 1, ..., n_u$

a) Compute the response $x_k^{iter}(t, \theta_k^{iter})$ to $u_k(t)$:

$$x_k^{iter}(t, \hat{\theta}_k^{iter}) = y^*(t) - \sum_{\substack{n=1 \\ n \neq k}}^{n_u} y_{u_n}(t, \hat{\theta}_n^{iter}) \quad (28)$$

b) Filter the noisy response of $x_k^{iter}(t)$, the auxiliary model outputs $y_{u_k}$ and the input signal $u_k$:

$$x_{k,f}^{iter}(t) = F_{svf}(p)x_k^{iter}(t) \quad (29)$$

$$y_{u_k,f}^{iter}(t) = F_{svf}(p)y_{u_k}^{iter}(t) \quad (30)$$

$$u_{k,f}(t) = F_{svf}(p)u_k(t) \quad (31)$$

c) Form the instrumental variable $\zeta_{k,f}$ and the regression $\varphi_{k,f}$ vectors:

$$\zeta_{k,f}(t) = \begin{bmatrix} p^{\beta_0} u_{k,f}(t) \dots p^{\beta_{M_k}} u_{k,f}(t) \\ -p^{\alpha_1} y_{u_k,f}^{iter}(t) \dots -p^{\alpha_{N_k}} y_{u_k,f}^{iter}(t) \end{bmatrix}^T \quad (32)$$

$$\varphi_{k,f}(t) = \begin{bmatrix} p^{\beta_0} u_{k,f}(t) \dots p^{\beta_{M_k}} u_{k,f}(t) \\ -p^{\alpha_1} x_{k,f}^{iter}(t) \dots -p^{\alpha_{N_k}} x_{k,f}^{iter}(t) \end{bmatrix}^T \quad (33)$$

d) Compute the IV-based solution $\hat{\theta}_k^{iter+1}$

$$\hat{\theta}_k^{iter+1} = [Z_{k,f} \Phi_{k,f}^T]^{-1} Z_{k,f} X_{k,f}^{iter} \quad (34)$$

with

$$Z_{k,f} = [\zeta_{k,f}(t_1), \dots, \zeta_{k,f}(t_H)] \quad (35)$$

$$\Phi_{k,f} = [\varphi_{k,f}(t_1), \dots, \varphi_{k,f}(t_H)] \quad (36)$$

$$X_{k,f}^{iter} = [x_{k,f}^{iter}(t_1), \dots, x_{k,f}^{iter}(t_H)] \quad (37)$$

e) Generate the auxiliary model output $y_{u_k}(t)$ using $\hat{\theta}_k^{iter+1}$:

$$y_{u_k}(t, \hat{\theta}_k^{iter+1}) = \frac{B_k(p, \hat{\theta}_k^{iter+1})}{A_k(p, \hat{\theta}_k^{iter+1})} u_k(t) \quad (38)$$

B) Iterate step (A) until:

$$test_k = \sum_{j=1}^{\dim(\hat{\theta}_k)} \left( \frac{\hat{\theta}_{k,j}^{iter+1} - \hat{\theta}_{k,j}^{iter}}{\hat{\theta}_{k,j}^{iter+1}} \right)^2 < \ell \quad (39)$$

where $\ell$ is a fixed edge, sufficiently small.

### Step3:

Compute the parametric covariance matrix $P_k$, after convergence.

for $k = 1, ..., n_u$

$$\hat{P}_k = \hat{\sigma}^2 \mathcal{H}_k^{-1} \quad (40)$$

with $\hat{\sigma}^2$ is the empirical estimation of the noise variance, given by:

$$\hat{\sigma}^2 = y^*(t) - \sum_{k=1}^{n_u} y_{u_k}^{iter}(t) \quad (41)$$

where $y_{u_k}^{iter}(t)$ is calculated at the last iteration.

$\mathcal{H}_k$ is the approximate Hessian given by:

$$\mathcal{H}_k = \frac{\partial \varepsilon^T}{\partial \theta_k} \frac{\partial \varepsilon}{\partial \theta_k} \quad (42)$$

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

$$\frac{\partial \varepsilon}{\partial \theta_k} = -\frac{\partial \hat{y}_{u_k}}{\partial \theta_k} = -\left[ \frac{\partial \hat{y}_{u_k}}{\partial b_{0,k}}, \dots, \frac{\partial \hat{y}_{u_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]^T \quad (43)$$

here $\frac{\partial \hat{y}_{u_k}}{\partial \hat{b}_{j,k}}$ and $\frac{\partial \hat{y}_{u_k}}{\partial a_{i,k}}$ are the output sensitivity functions, given by:

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

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

An algorithm for better set the cut-off frequency $\omega_c$, on the basis of the Shannon's frequency $\omega_s$, of the filter has been implemented (equation 19).

## VI. NUMERICAL EXAMPLE

Consider the following MISO system. Data are generated using the following system:

$$S: \begin{cases} y_u(t) = \frac{1}{3p^{0.5+1}} u_1(t) + \frac{2}{2p^{0.5+1}} u_2(t) + \frac{5}{p^{0.5+1}} u_3(t) \\ y(t_h) = y_u(t_h) + e(t_h) \end{cases} \quad (46)$$

Input signals $u_i$, $i = 1, 2, 3$ are three uncorrelated pseudo-random-binary-sequences (prbs). The prbs amplitudes are fixed between -1 and 1. The prbs sequences are set such as $3 \times \tau_{max}$ (with: $\tau_{max} = \max(\tau_1, \tau_2, \tau_3)$); where $\tau_{i=1,2,3}$ are the time constants of the system 46. The sampling period is chosen as $T = 0.05s$ in those simulations.
## A. White noise

The noise-free output $y_u(t)$, which is the sum of output responses is corrupted by a zero mean Gaussian white noise. Two signal to noise ratio (SNR) has been tested to validate the two methods $SNR_1 = 30dB$ and $SNR_2 = 10dB$.

All the differential orders are supposed known and set to $\nu_1 = \nu_2 = \nu_3 = 0.5$. The unknown coefficients are estimated by applying the **oe** and **ivsvf** methods. A Monte Carlo analysis, for $N_{exp} = 200$ experiments, is simulated in this section to demonstrate the efficiency of the MISO version of **ivsvf** and **oe** methods.

The distributions of estimates obtained by applying the **oe** and **ivsvf** methods are shown respectively in Fig.2 and Fig.3. Table I and Table II summarize the means of estimates and their standard deviations.

For both **oe** and **ivsvf**, the estimates are unbiased and their distributions are close to the true parameters.

### 1) Additive noise (SNR = 30dB):

<table><caption>TABLE I
MONTE CARLO SIMULATION WITH 200 RUNS ($\hat{\theta}$ BEING THE MEAN, $\sigma(\hat{\theta})$ BEING THE STANDARD DEVIATION OF THE PARAMETER ESTIMATES) FOR WHITE NOISE SNR = 30dB.</caption><thead><tr><th rowspan="2">True</th><th colspan="2">ivsvf</th><th colspan="2">oe</th></tr><tr><th>θ̂</th><th>σ(θ̂)</th><th>θ̂</th><th>σ(θ̂)</th></tr></thead><tbody><tr><td>b<sub>0,1</sub> 1</td><td>1.0028</td><td>0.0008</td><td>1.0064</td><td>0.0008</td></tr><tr><td>a<sub>1,1</sub> 3</td><td>3.0311</td><td>0.0126</td><td>3.0324</td><td>0.0131</td></tr><tr><td>b<sub>0,2</sub> 2</td><td>1.9996</td><td>0.0004</td><td>2.0006</td><td>0.0004</td></tr><tr><td>a<sub>1,2</sub> 2</td><td>1.9997</td><td>0.0008</td><td>1.9994</td><td>0.0008</td></tr><tr><td>b<sub>0,3</sub> 5</td><td>5.0137</td><td>0.0001</td><td>5.0044</td><td>0.0001</td></tr><tr><td>a<sub>1,3</sub> 1</td><td>1.0127</td><td>0.0001</td><td>1.0012</td><td>0.0001</td></tr></tbody></table>

<table><caption>TABLE II
MONTE CARLO SIMULATION WITH 200 RUNS ($\hat{\theta}$ BEING THE MEAN, $\sigma(\hat{\theta})$ BEING THE STANDARD DEVIATION OF THE PARAMETER ESTIMATES) FOR WHITE NOISE SNR = 10dB.</caption><thead><tr><th rowspan="2">True</th><th colspan="2">ivsvf</th><th colspan="2">oe</th></tr><tr><th>θ̂</th><th>σ(θ̂)</th><th>θ̂</th><th>σ(θ̂)</th></tr></thead><tbody><tr><td>b<sub>0,1</sub> 1</td><td>1.1158</td><td>0.0388</td><td>1.1247</td><td>0.0520</td></tr><tr><td>a<sub>1,1</sub> 3</td><td>3.5208</td><td>0.6294</td><td>3.5283</td><td>0.7152</td></tr><tr><td>b<sub>0,2</sub> 2</td><td>2.0023</td><td>0.0294</td><td>2.0034</td><td>0.0433</td></tr><tr><td>a<sub>1,2</sub> 2</td><td>2.0054</td><td>0.0610</td><td>2.0060</td><td>0.0857</td></tr><tr><td>b<sub>0,3</sub> 5</td><td>5.0191</td><td>0.0079</td><td>5.0180</td><td>0.0079</td></tr><tr><td>a<sub>1,3</sub> 1</td><td>1.0137</td><td>0.0009</td><td>1.0102</td><td>0.0009</td></tr></tbody></table>

### 2) Additive noise (SNR = 10dB):

<figure><img src="image_3.png" /><figcaption>Fig. 2. The set of the Monte Carlo analysis with white additive noise (SNR = 30dB).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. The set of the Monte Carlo analysis: white additive noise (SNR = 10dB).</figcaption></figure>

## B. Colored noise

In this section, the measurements of output signal are supposed contaminated by a colored noise. It is generated by a low-pass filter excited by a white noise with $SNR = 20dB$.

$N_{exp} = 200$ runs of different colored noises are generated to analyse the robustness of the developed methods to additive noises. Both **oe** and **ivsvf** methods are applied to estimate the coefficients of the fractional MISO model described by (45).

The estimates are unbiased and their distributions are
<figure><img src="image_5.png" /><figcaption>Fig. 4. The set of the Monte Carlo analysis with colored noise ($SNR = 20dB$)</figcaption></figure>

centered on the true values (see Fig.4 and Table III).

<table><caption>TABLE III
MONTE CARLO SIMULATION WITH 200 RUNS ($\hat{\theta}$ BEING THE MEAN, $\sigma(\hat{\theta})$ BEING THE STANDARD DEVIATION OF THE PARAMETER ESTIMATES) FOR COLORED NOISE $SNR = 20dB$.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">ivsvf</th><th colspan="2">oe</th></tr><tr><th>θ̂</th><th>σ(θ̂)</th><th>θ̂</th><th>σ(θ̂)</th></tr></thead><tbody><tr><td>b<sub>0,1</sub></td><td>1</td><td>1.0030</td><td>0.0080</td><td>1.0049</td><td>0.0080</td></tr><tr><td>a<sub>1,1</sub></td><td>3</td><td>2.9940</td><td>0.2056</td><td>3.0084</td><td>0.2121</td></tr><tr><td>b<sub>0,2</sub></td><td>2</td><td>1.9983</td><td>0.0040</td><td>2.0023</td><td>0.0044</td></tr><tr><td>a<sub>1,2</sub></td><td>2</td><td>2.0003</td><td>0.0147</td><td>2.0018</td><td>0.0150</td></tr><tr><td>b<sub>0,3</sub></td><td>5</td><td>5.0151</td><td>0.0067</td><td>5.0035</td><td>0.0048</td></tr><tr><td>a<sub>1,3</sub></td><td>1</td><td>1.0125</td><td>0.0009</td><td>1.0011</td><td>0.0009</td></tr></tbody></table>

## VII. CONCLUSION

In this paper, two methods have been extended to the fractional multiple-input-single-output MISO models, the first one is the instrumental variable state-variable filter ivsvf approach and the second one is the output-error oe approach. A comparison between the two methods has been introduced in this paper. A Monte carlo simulations where used to validate the two new methods. As it is illustrated in the numerical results the ivsvf method shows a better estimation than the oe method despite of the noise ratio presence.

## REFERENCES

[1] P.C. Young. An instrumental variable method for real-time identification of a noisy process. Automatica, 6:271-287, 1970

[2] P.C. Young and A.J. Jakeman. Refined instrumental variable methods of time-series analysis: Part iii, extensions. International Journal of Control, 31:741-764, 1980.

[3] P.C. Young. Recursive estimation and Time-Series Analysis. Springer-Verlag, Berlin, 2nd edition, 2011.

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

[5] H. Garnier, M. Gilson, P.C. Young, and E. Huselstein. An optimal iv technique for identifying continuous-time transfer function model of multiple input systems. Control Engineering Practice, 15(4):471-486, 2007.

[6] R. S. Blom and P. M. J. Van den Hof. Multivariable frequency domain identification using IV-based linear regression. 49th IEEE Conference on Decision and Control (CDC), 1148-1153, 2010.

[7] O. Cois, A. Oustaloup, E. Battaglia, and J.-L. Battaglia. Non integer model from modal decomposition for timedomain system identification. In 12th IFAC Symposium on System Identification, SYSID, Santa Barbara, USA, June 2000.

[8] M. Aoun, R. Malti, F. Lebron, and A. Oustaloup. Numerical simulations of fractional systems. In ASMEIDETC/CIE Conferences, pages DETC2003/VIB-48389, Chicago Il USA, 9 2003.

[9] R. Malti, S. Victor, A. Oustaloup, and H. Garnier. An optimal instrumental variable method for continuous-time fractional model identification. In 17th IFAC World Congress, pages 14379-14384, Seoul South Korea, 2008.

[10] O. Cois, A. Oustaloup, T. Poinot, and J.-L. Battaglia. Fractional state variable filter for system identification by fractional model. In 6th IEEE ECC, Porto, Portugal, Sep. 2001.

[11] M. Chetoui, R. Malti, M. Thomassin, M. Aoun, S. Najar, and M. N. Abdelkrim. Third-order cumulants based method for continuous-time errors-in-variables system identification by fractional models. In Eighth International Multi-Conference on Systems, Signals Devices (pp. 1-6). IEEE. 2011.

[12] M. Chetoui, M. Thomassin, R. Malti, M. Aoun, S. Najar, and M. N. 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(5), 860-872. 2013.

[13] Z. Yakoub, M. Amairi, M. Chetoui and M. Aoun. On the closed-loop system identification with fractional models. Circuits, Systems, and Signal Processing, 34(12), 3833-3860. 2015.

[14] M. Aoun, R. Malti, F. Lebron and A. Oustaloup. Synthesis of fractional Laguerre basis for system approximation. Automatica, 43(9), 1640-1648. 2007.

[15] R. Malti, M. Aoun, J. Sabatier and A. Oustaloup. Tutorial on system identification using fractional differentiation models. IFAC Proceedings Volumes, 39(1), 606-611. 2006.

[16] T. Salem, M. Chetoui, and M. Aoun. Instrumental variable based methods for continuous-time linear parameter varying system identification with fractional models. In 2016 24th Mediterranean Conference on Control and Automation (MED) (pp. 640-645). IEEE. 2016.

[17] F. Chen, M. Gilson, H. Garnier and T. Liu. Robust time-domain output error method for identifying continuous-time systems with time delay. Systems Control Letters, 102, 81-92. 2017. doi:10.1016/j.sysconle.2017.01.009.