# Instrumental Variable Identification of Hybrid Fractional Box-Jenkins Models

Stéphane Victor, Rachid Malti, Pierre Melchior, Alain Oustaloup

► To cite this version:

Stéphane Victor, Rachid Malti, Pierre Melchior, Alain Oustaloup. Instrumental Variable Identification of Hybrid Fractional Box-Jenkins Models. 18th IFAC World Congress, Aug 2011, Milan, Italy. pp.4314-4319, 10.3182/20110828-6-IT-1002.01107. [https://hal-05454216](https://hal.science/hal-05454216v1)

HAL Id: hal-05454216

Submitted on 12 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Instrumental Variable Identification of Hybrid Fractional Box-Jenkins Models

Stéphane Victor, Rachid Malti, Pierre Melchior, and Alain Oustaloup

Université de Bordeaux - IPB. IMS, UMR 5218 CNRS, 351 cours de la Libération, F 33405 Talence cedex, France - {stephane.victor, rachid.malti, pierre.melchior, alain.oustaloup}@laps.ims-bordeaux.fr

**Abstract:** this paper deals with continuous-time system identification using fractional differentiation models in a colored noisy output context. An optimal instrumental variable method for identifying fractional models in a hybrid Box-Jenkins form is described. The relationship between the measured input and the output is a fractional continuous-time transfer function, and the noise is a discrete-time AR or ARMA process.

**Keywords:** fractional continuous-time models; fractional differentiation; instrumental variable; state variable filter; system identification; time-domain.

## 1. INTRODUCTION

Instrumental Variable (iv) approaches to system identification and parameter estimation have been developed since 1970s as detailed in the survey paper by Young and Jakeman (1980) and all references therein.

The Simplified Refined Instrumental Variable (srivc) algorithm for fractional systems, namely srivcf, has been developed in Malti et al. (2008). In that paper, the least-square (ls) method and the instrumental variable (iv) combined with state variable filter (svf) are recalled. The model structure is fixed and only the coefficients were identified in a noisy context. An iterative optimal instrumental variable method for fractional continuous-time model has been implemented where the auxiliary model and the state variable filters are updated at each iteration. This optimal iv method presents the advantage of yielding asymptotically consistent estimates in the presence of additive noise and of minimum variance in the presence of white noise. However, if the noise is not white, the variance is no more minimum.

In this paper, the dimension of the true system and the commensurate<fnref n="1" /> order are assumed to be known. Although the iv aspects ensure that the parameter estimates are asymptotically unbiased, even in noisy colored context, the estimates are not statistically efficient in this situation because the svf are not suited in colored noise context. Recently, Refined Instrumental Variables for Continuous-time system (Young et al. (2006)) has been developed for identifying hybrid continuous-time Box-Jenkins models, where the system model is continuous-time (CT) and the noise model is discrete-time (DT). Here, the CT Box and Jenkins models are extended to fractional models: a fractional CT model links the input to the output and a DT ARMA process is used for noise modeling.

A fractional mathematical model is based on a fractional differential equation:

$$y(t) + a_1 p^{\alpha_1} y(t) + \dots + a_{m_A} p^{\alpha_{m_A}} y(t) = b_0 p^{\beta_0} u(t) + b_1 p^{\beta_1} u(t) + \dots + b_{m_B} p^{\beta_{m_B}} u(t), \quad (1)$$

where $(a_j, b_i) \in \mathbb{R}^2$, $p = \frac{d}{dt}$ is the differential operator, the differentiation orders $0 < \alpha_1 < \alpha_2 < \dots < \alpha_{m_A}$, $0 < \beta_0 < \beta_1 < \dots < \beta_{m_B}$ are allowed to be non-integer positive numbers. The concept of differentiation to an arbitrary order (non-integer) $p^\gamma$, with $\gamma \in \mathbb{R}_+^*$ (set of strictly positive real numbers) was defined in the 19th century by Riemann and Liouville. The $\gamma$-order fractional derivative of $y(t)$ is defined as being an integer derivative of order $\lfloor \gamma \rfloor + 1$ ($\lfloor \cdot \rfloor$ stands for the floor operator) of a non-integer integral of order $\gamma - \lfloor \gamma \rfloor$ Samko et al. (1993):

$$p^{\gamma} y(t) = p^{\lfloor \gamma \rfloor + 1} \left( \mathbf{I}^{\lfloor \gamma \rfloor + 1 - \gamma} y(t) \right) \triangleq p^{\lfloor \gamma \rfloor + 1} \left( \frac{1}{\Gamma(\lfloor \gamma \rfloor + 1 - \gamma)} \int_0^t \frac{y(\tau) d\tau}{(t - \tau)^{\gamma - \lfloor \gamma \rfloor}} \right), \quad (2)$$

where I stands for the integral operator and the Euler $\Gamma$ function is defined in the set of real numbers except zero and negative whole numbers ($x \in \mathbb{R}^* \setminus \mathbb{N}^-$) as:

$$\Gamma(x) = \int_{0}^{\infty} e^{-t} t^{x-1} dt. \quad (3)$$

The Laplace transform is a more concise algebraic tool generally used to represent fractional systems, Oldham and Spanier (1974):

$$\mathcal{L}\{y^{(\gamma)}(t)\} = s^{\gamma} X(s) \quad \text{if } y(t) = 0 \ \forall t \le 0. \quad (4)$$

This property allows to write the fractional differential equation (1) in a transfer function form:

$$G(s) = \frac{B(s)}{A(s)} = \frac{\sum_{i=0}^{m_B} b_i s^{\beta_i}}{1 + \sum_{j=1}^{m_A} a_j s^{\alpha_j}}. \quad (5)$$

<footnote n="1">¹ All differentiation orders are exactly divisible by the same number, an integral number of times, and the biggest number is always chosen.</footnote>
**Stability theorem (Matignon (1998), extended):** A commensurate transfer function $G(s) = \frac{T(s^\gamma)}{R(s^\gamma)}$, T and R being coprime polynomials, is BIBO stable if and only if:

$$0 < \gamma < 2 \qquad (6)$$

and

$$|\arg(s_k)| > \gamma \frac{\pi}{2}, \quad \forall s_k \in \mathbb{C} / R(s_k) = 0. \qquad (7)$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Stability region: a system is stable iff its $s^\gamma$-poles are inside the green area ($\gamma < 1$)</figcaption></figure>

The aim is to develop a method for time-domain system identification workable with any fractional-model-time-domain simulation algorithm. Here, the recursive approximation of a fractional differentiation with poles and zeros developed by Oustaloup (1995) is used for simulating the fractional operator.

## 2. SYSTEM IDENTIFICATION

Most of the direct CT model identification techniques do not consider the estimation of a noise model leading into statistically inefficient estimates (high estimate variance).

The **Refined Instrumental Variables** for **Continuous system (rivc**: see Young and Jakeman (1980)) method and its **Simplified version srivc** (Young (2002)), when the additive measurement noise is assumed to be white, are the only **iv** methods that can be interpreted in optimal statistical terms. A logical extension of srivc for fractional models, srivcf, has been developed in Malti et al. (2008).

This approach involves an iteratively-adaptive prefiltering method based on a quasi-optimal statistical solution to the problem in this white noise case. This technique has the advantage that the prefilters are automatically optimized in an iterative algorithm. Optimal instrumental variables are obtained from an iteratively updated auxiliary model. At the first iteration, the auxiliary model is initialized using a least square estimate obtained with state variable filters. The initial $\hat{A}_0(s)$ polynomial does not have to be particularly accurate provided the prefilter is stable and does not seriously attenuate input/output signals within the system bandwidth.

### 2.1 Problem formulation

Consider $K + 1$ regular samples of input $u(t)$ and noisy output $y^*(t)$ data collected from $t_0 = 0$ to $t_K = KT_s$.

In order to have good statistical estimation results, it is necessary to consider stochastic errors that affect the

output signal. When the additive measurement noise $\xi(t)$ is colored, the model describing the system is of Box-Jenkins type:

$$\left\{ \begin{aligned} y(t) &= G(p)u(t), \\ \xi(t) &= H(p)e(t), \\ y^*(t) &= y(t) + \xi(t), \end{aligned} \right. \qquad (8)$$

where $y$ is the hypothetical noise-free deterministic system output, the input/output model $G(p)$ is described by (5) and the noise model $H(p)$, supposed stable and invertible, is described by an Auto-Regressive (AR) or Auto-Regressive Moving Average (ARMA) process:

$$H(p) = \frac{C(p)}{D(p)}, \qquad (9)$$

and where $e$ is a zero-mean white additive gaussian noise.

Since the data are sampled and the representation of a white noise is in discrete-time, a more suited noise model is a discrete-time ARMA($r, v$) process:

$$\mathcal{H}(q^{-1}) = \frac{C(q^{-1})}{D(q^{-1})} = \frac{1 + \sum_{i=0}^{v} c_i q^{-i}}{1 + \sum_{j=1}^{r} d_j q^{-j}}, \qquad (10)$$

where $q^{-1}$ is the backward shift operator such that $q^{-1}y(t_k) = y(t_{k-1})$.

A hybrid model is then formulated with

$$\left\{ \begin{aligned} y(t) &= G(p)u(t) \\ \xi(t_k) &= \mathcal{H}(q^{-1})e(t_k) \\ y^*(t_k) &= y(t_k) + \xi(t_k), \end{aligned} \right. \qquad (11)$$

where $e(t_k)$ is the white gaussian noise sampled at $t_k$, the system model is defined in continuous-time as in (5) and the noise model is discrete-time as in (10).

When the differentiation orders are known, the parameter vector

$$\theta = [\rho, \eta]^T, \qquad (12)$$

is composed of the $m_A + m_B + 1$ coefficients of the system model

$$\rho = [b_0, b_1, \dots, b_{m_B}, a_1, \dots, a_{m_A}]^T,$$

and of the $r + v + 1$ coefficients of the noise model

$$\eta = [c_0, c_1, \dots, c_v, d_1, \dots, d_r]^T.$$

### 2.2 Least square estimate

In white gaussian additive noise ($H(p) = 1$), the measured output is corrupted by an additive measurement noise $e$, $y^*(t) = y(t) + e(t)$, and the equation error can be formulated as:

$$\varepsilon(t) = y^*(t) - \varphi^*(t)^T \theta, \qquad (13)$$

where

$$\varphi^*(t)^T = \begin{bmatrix} u^{(\beta_0)}(t) \cdots u^{(\beta_{m_B})}(t) \\ -y^{*(\alpha_1)}(t) \cdots -y^{*(\alpha_{m_A})}(t) \end{bmatrix} \qquad (14)$$

and

$$\theta^T = \rho^T = [b_0 \ b_1 \ \dots \ b_{m_B} \ a_1 \ \dots \ a_{m_A}]. \qquad (15)$$

Minimizing the $L_2$ norm of $\varepsilon(t)$ with respect to $\theta$, leads to the following least square estimates:

$$\hat{\theta}_{ls} = \left[ \int_{t_0}^{t_K} \varphi^*(t)^T \varphi^*(t) dt \right]^{-1} \int_{0}^{t_K} \varphi^*(t)^T y^*(t) dt. \qquad (16)$$
After a discretization, by defining $\mathbf{Y}$ as a column vector of the system output and $\mathbf{\Phi}^*$ as a regression matrix:

$$ \mathbf{\Phi}^* = [\varphi^*(t_0) \quad \varphi^*(t_1) \quad \dots \quad \varphi^*(t_K)]^T, \quad (17) $$

$\hat{\theta}_{ls}$ can be approximated by:

$$ \hat{\theta}_{ls} = (\mathbf{\Phi}^{*\mathrm{T}}\mathbf{\Phi}^{*})^{-1} \mathbf{\Phi}^{*\mathrm{T}}\mathbf{Y}^{*}. \quad (18) $$

Direct fractional differentiations of noisy output lead however to inaccurate results. Thus, *srivcf* is preferred to the direct *ls* estimate.

## 2.3 Simplified refined instrumental variable for continuous-time fractional models (srivcf)

The *srivc* stochastic approach proposed by Young (2002) includes an adaptive prefiltering based on the optimal solution when the noise $\xi$ is a white gaussian one ($H(p) = 1$ and the parameter vector is resumed to $\theta = \rho$). This method is a logical extension of the instrumental variable estimator with state variable filters (Cois et al. (2001)). This estimator extended to fractional models as *srivcf* (*srivc* for fractional models) is recalled (see Malti et al. (2008)).

Following the Prediction Error Minimization (PEM) approach, a suitable error function $\varepsilon(t)$ in CT is given by the output error (OE):

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

The minimization of the quadratic criterion function $\varepsilon(t)$ provides the basis for stochastic estimation. The error function can be rewritten as:

$$ \begin{aligned} \varepsilon(t) &= \frac{1}{A(p)} [A(p) y^*(t) - B(p) u(t)] \\ &= A(p) \left( \frac{1}{A(p)} y^*(t) \right) - B(p) \left( \frac{1}{A(p)} u(t) \right). \end{aligned} \quad (20) $$

By considering $F(p) = \frac{1}{A(p)}$ as the filter, $y_f(t) = F(p)y^*(t)$ and $u_f(t) = F(p)u(t)$, the error becomes

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

or,

$$ \begin{aligned} \varepsilon(t) = & y_f^*(t) + a_1 y_f^{*(\alpha_1)}(t) + \dots + a_{m_A} y_f^{*(\alpha_{m_A})}(t) \\ & - b_0 u_f^{(\beta_0)}(t) - b_1 u_f^{(\beta_1)}(t) - \dots - b_{m_B} u_f^{(\beta_{m_B})}(t), \end{aligned} \quad (22) $$

where

$$ \left\{ \begin{aligned} u_f^{\beta_i}(t) &= p^{\beta_i} F(p) u(t), & i &= 0, \dots, m_B \\ y_f^{\alpha_j}(t) &= p^{\alpha_j} F(p) y^*(t), & j &= 1, \dots, m_A. \end{aligned} \right. \quad (23) $$

In practical cases, $A(p)$ is not known. The filter $F(p)$ is then updated by an iterative procedure which estimates the denominator $\hat{A}(p)$. This algorithm enables to optimize the instruments. Initialized by the $\hat{\theta}_{iv/svf}$ estimates (see Malti et al. (2008) for the details), as the first estimation of the denominator is available, the new derivatives of input and output signals are computed with the set of filters:

$$ F_{\mathrm{iter}}(p) = \frac{1}{\hat{A}(p)} = \frac{1}{1 + \sum_{j=1}^{m_A} \hat{a}_j p^{\alpha_j}}, \quad (24) $$

where iter = 1, 2, ... stands for the iteration number and $\hat{a}_j$ the estimated $a_j$ coefficient at the iteration iter.

Instruments and input/output filtered derivatives are computed at each iteration by:

$$ \left\{ \begin{aligned} u_f^{(\beta_i)}(t) &= p^{\beta_i} F_{\mathrm{iter}}(p) u(t) \\ y_f^{*(\alpha_j)}(t) &= p^{\alpha_j} F_{\mathrm{iter}}(p) y^*(t) \\ y_f^{iv(\alpha_j)}(t) &= p^{\alpha_j} F_{\mathrm{iter}}(p) y^{iv}(t), \end{aligned} \right. \quad (25) $$

where $y^{iv}(t) = \frac{\hat{B}(p)}{A(p)}u(t)$. The regression vectors $\varphi_f^*(t)$ and $\varphi_f^{iv}(t)$ are formed as in (14):

$$ \begin{aligned} \varphi_f^*(t)^T &= \begin{bmatrix} u_f^{(\beta_0)}(t) \cdots u_f^{(\beta_{m_B})}(t) \\ -y_f^{*(\alpha_1)}(t) \cdots -y_f^{*(\alpha_{m_A})}(t) \end{bmatrix} \\ \varphi_f^{iv}(t)^T &= \begin{bmatrix} u_f^{\beta_0}(t) \cdots u_f^{(\beta_{m_B})}(t) \\ -y_f^{iv(\alpha_1)}(t) \cdots -y_f^{iv(\alpha_{m_A})}(t) \end{bmatrix}. \end{aligned} \quad (26) $$

Finally, the parameter vector $\hat{\theta}_{\mathrm{iter}}^{srivcf}$ is computed as:

$$ \hat{\theta}_{\mathrm{iter}}^{srivcf} = (\mathbf{\Phi}_{f}^{iv\mathrm{T}}\mathbf{\Phi}_{f}^{*})^{-1} \mathbf{\Phi}_{f}^{iv\mathrm{T}}\mathbf{Y}_{f}^{*}, \quad (27) $$

with

$$ \mathbf{Y}_f^* = [y_f^*(t_0) \quad y_f^*(t_1) \quad \dots \quad y_f^*(t_K)]^T, $$

$$ \mathbf{\Phi}_f^* = [\varphi_f^*(t_0) \quad \varphi_f^*(t_1) \quad \dots \quad \varphi_f^*(t_K)]^T, $$

$$ \mathbf{\Phi}_f^{iv} = [\varphi_f^{iv}(t_0) \quad \varphi_f^{iv}(t_1) \quad \dots \quad \varphi_f^{iv}(t_K)]^T. \quad (28) $$

As in the rational case, the *srivcf* estimator is asymptotically unbiased whatever additive zero-mean noise is present. When the additive noise is white, *srivcf* has minimum variance and the asymptotic covariance matrix of the estimation errors associated to the estimate of $\hat{\theta}^{srivcf}$:

$$ \mathbf{P}_{\hat{\theta}} = \hat{\sigma}^2 (\mathbf{\Phi}_f^{iv\mathrm{T}} \mathbf{\Phi}_f^{*})^{-1}, \quad (29) $$

where $\hat{\sigma}^2$ is the empirical estimation of noise variance.

The associated estimation model can be written as

$$ y_f(t) = \varphi_f(t)^T \theta + \varepsilon(t_k), \quad (30) $$

where $\varepsilon(t_k)$ is the equation error, and

$$ \varphi_f(t)^T = \begin{bmatrix} u_f^{(\beta_0)} u_f(t) \cdots u_f^{(\beta_{m_B})}(t) \\ -y_f^{(\alpha_1)}(t) \cdots -y_f^{(\alpha_{m_A})}(t) \end{bmatrix}. \quad (31) $$

Hence, the estimation model (30) forms a basis for the definition of a likelihood function and maximum likelihood estimation.

## Summary of the srivcf algorithm

**Step 1** Use the *iv* method combined to fractional state variable filter to generate an initial (iter = 1) estimate of the transfer function model parameter vector $\hat{\theta}^{(1)}$.

**Step 2** Iterative *iv* estimation with prefilters.

for iter = 2 to convergence

(i) Generate the *iv* $y^{iv}$ from the auxiliary model with the estimated polynomials based on the estimated parameter vector $\hat{\theta}_{\mathrm{iter}-1}^{srivcf}$.

(ii) Update the filter $F_{\mathrm{iter}}(s, \hat{\theta})$ in (24) with the new estimated parameters. Then, evaluate the prefiltered derivatives of $u(t)$, $y(t)$ and $y^{iv}(t)$ as in (25).

(iii) Based on these prefiltered data, compute the new estimates $\hat{\theta}_{\mathrm{iter}}^{srivcf}$ as in (27).

end for
**Step 3** Compute the estimated parametric error covariance matrix with (29).

It may happen, when the bias is important, that the estimated parameters from Step 1 lead to unstable models. In this case, an empirical method<fnref n="2" /> can be used to stabilize unstable poles.

## 2.4 Refined instrumental variable for continuous-time fractional Box-Jenkins models (rivcf)

If the noise $\xi(t)$ is colored ($H(p) \neq 1$ and the parameter vector $\theta$ is defined by $[\rho, \eta]^T$), then the *srivcf* method has no claim to optimality in statistical terms as the variance of the estimated parameters is not minimal. However, thanks to the *iv* mechanism, the estimates are consistent (asymptotically unbiased). The algorithm is very robust and normally produces estimates with reasonable statistical efficiency, meaning that the variance is low but not minimum. The algorithm is extended to the colored noise case, namely *rivcf* (rivc for fractional systems, see Garnier and Wang (2008) for the rational case) method, by using a hybrid approach: the system model is a fractional CT transfer function and the noise model is a DT AR or ARMA process.

By following the usual PEM approach in the hybrid case, which corresponds to the maximum likelihood estimation given the gaussian hypothesis on $e$, the error function at $t_k$ is given by

$$e(t_k) = \frac{\mathcal{D}(q^{-1})}{\mathcal{C}(q^{-1})} \left\{ y^*(t_k) - \frac{B(p)}{A(p)} u(t_k) \right\}, \quad (32)$$

which can also be written as

$$e(t_k) = \frac{\mathcal{D}(q^{-1})}{\mathcal{C}(q^{-1})} \left\{ \frac{1}{A(p)} [A(p)y^*(t_k) - B(p)u(t_k)] \right\}, \quad (33)$$

where the DT prefilter $\frac{\mathcal{D}(q^{-1})}{\mathcal{C}(q^{-1})}$ corresponds to the inverse of the noise model ARMA($r, v$).

In the last two equations, the DT and CT operators are used in order to indicate the hybrid nature of the estimation problem. Thus, the operations such that $\frac{B(p)}{A(p)}u(t_k)$ imply that the variable $u(t_k)$ needs to be interpolated.

The *ls* error criterion minimisation on $e(t_k)$ being at the heart of the stochastic estimation, the equation (33) can be considered under the alternative form

$$e(t_k) = A(p)y_f^*(t_k) - B(p)u_f(t_k), \quad (34)$$

where $y_f^*(t_k)$ and $u_f(t_k)$ are the sampled outputs of the hybrid prefiltering operations, namely a CT filtering of $y^*$ and $u$ with the system model denominator:

$$F_c^{\text{opt}}(p) = \frac{1}{A(p)}, \quad (35)$$

and a DT filtering of the error $e$ with the inverse noise model:

$$F_d^{\text{opt}}(q^{-1}) = \frac{\mathcal{D}(q^{-1})}{\mathcal{C}(q^{-1})}. \quad (36)$$

The system and noise models being unknown in practical cases, iterative procedure are suggested in the *iv* estimation methods. The optimal prefilters in (35) and in (36) are respectively replaced by

$$F_{c, \text{iter}}(p) = \frac{1}{\hat{A}(p)} \quad (37)$$

and

$$F_{d, \text{iter}}(q^{-1}) = \frac{\hat{D}(q^{-1})}{\hat{C}(q^{-1})}, \quad (38)$$

computed iteratively where $\hat{A}$, $\hat{C}$ and $\hat{D}$ correspond respectively to the estimation of $A$ obtained at iteration (iter - 1), and to the estimations of $C$ and $D$ computed at the current iteration (iter). Moreover, the instrumental variable is generated from the CT auxiliary model

$$y^{\text{iv}}(t) = \frac{\hat{A}(p)}{\hat{B}(p)} u(t). \quad (39)$$

The filtered derivatives of the input, the output and the instrumental variable are computed at each iteration by:

$$\left\{ \begin{aligned} u_f^{(\beta_i)}(t_k) &= F_{d, \text{iter}}(q^{-1}) u_{fc}^{(\beta_i)}(t_k), & i = 0, \dots, m_B, \\ y_f^*(\alpha_j)(t_k) &= F_{d, \text{iter}}(q^{-1}) y_{fc}^*(\alpha_j)(t_k), & j = 1, \dots, m_A, \\ y_f^{\text{iv}}(\alpha_j)(t_k) &= F_{d, \text{iter}}(q^{-1}) y_{fc}^{\text{iv}}(\alpha_j)(t_k), & j = 1, \dots, m_A, \\ u_{fc}^{(\beta_i)}(t) &= p^{\beta_i} F_{c, \text{iter}}(p) u(t), & i = 0, \dots, m_B, \\ y_{fc}^*(\alpha_j)(t) &= p^{\alpha_j} F_{c, \text{iter}}(p) y^*(t), & j = 1, \dots, m_A, \\ y_f^{\text{iv}}(\alpha_j)(t) &= p^{\alpha_j} F_{c, \text{iter}}(p) y^{\text{iv}}(t), & j = 1, \dots, m_A. \end{aligned} \right.$$

The regression vector $\varphi_f^*(t_k)$ and the instrumental vector $\varphi_f^{\text{iv}}(t_k)$, which depend on $\hat{\rho}_{\text{iter}-1}$ and $\hat{\eta}_{\text{iter}}$, are computed:

$$\begin{aligned} \varphi_f^*(t)^T &= \begin{bmatrix} u_f^{(\beta_0)}(t) \cdots u_f^{(\beta_{m_B})}(t) \\ -y_f^*(\alpha_1)(t) \cdots -y_f^*(\alpha_{m_A}) (t) \end{bmatrix} \\ \varphi_f^{\text{iv}}(t)^T &= \begin{bmatrix} u_f^{\beta_0}(t) \cdots u_f^{(\beta_{m_B})}(t) \\ -y_f^{\text{iv}}(\alpha_1)(t) \cdots -y_f^{\text{iv}}(\alpha_{m_A}) (t) \end{bmatrix}. \end{aligned} \quad (40)$$

Finally, the parameter vector $\hat{\rho}_{\text{iter}}^{\text{rivcf}}$ is calculated as:

$$\hat{\rho}_{\text{iter}}^{\text{rivcf}} = \left[ \mathbf{\Phi}_f^{\text{iv}^T} \mathbf{\Phi}_f^* \right]^{-1} \mathbf{\Phi}_f^{\text{iv}^T} \mathbf{Y}_f^*, \quad (41)$$

with

$$\begin{aligned} \mathbf{Y}_f^* &= [y_f^*(t_0), y_f^*(t_1), \dots, y_f^*(t_K)]^T, \\ \mathbf{\Phi}_f^* &= [\varphi_f^*(t_0), \varphi_f^*(t_1), \dots, \varphi_f^*(t_K)]^T, \\ \mathbf{\Phi}_f^{\text{iv}} &= [\varphi_f^{\text{iv}}(t_0), \varphi_f^{\text{iv}}(t_1), \dots, \varphi_f^{\text{iv}}(t_K)]^T. \end{aligned} \quad (42)$$

The parameter vector $\hat{\eta}_{\text{iter}}$ is estimated at each iteration by applying any ARMA model estimation algorithm (Ljung (1999); Söderström and Stoica (1989)) on the output error, assimilated to the colored noise at each iteration:

$$\hat{\xi}(t_k) = y^*(t_k) - y^{\text{iv}}(t_k). \quad (43)$$

This algorithm, based on linear regression, is defined from a quadratic error criterion that is convex. Therefore, in all identifiable situations, the *rivcf* algorithm converges to the true parameters and the estimates will be consistent and asymptotically efficient.

Asymptotic covariance matrix estimation $\mathbf{P}_{\rho}$ associated to the estimates $\hat{\rho}_{\text{iter}}^{\text{rivcf}}$, obtained at the convergence of $\hat{\rho}_{\text{iter}}^{\text{rivcf}}$, is given by:

<footnote n="2">² This method consists of replacing every unstable $s^\gamma$-poles by its symmetrical one with respect to its closest stability axis (Fig. 1)</footnote>
$$ \mathbf{P}_{\hat{\rho}^{rivcf}} = \hat{\sigma}^2 \left[ \sum_{k=0}^{K} \varphi_f^{iv}(t_k) \varphi_f^{iv}(t_k)^{\mathrm{T}} \right]^{-1}, \quad (44) $$

where $\varphi_f^{iv}(t_k)$ is the instrumental vector at convergence and $\hat{\sigma}^2$ is the empirical estimation of noise variance.

## Summary of the full rivcf Algorithm

### Step 1 Initialization

Apply the *rivcf* algorithm in order to obtain an initial estimation at iteration iter = 1 of the parameter vector $\hat{\rho}_1$, this estimation being used to generate the CT initial prefilter $F_{c,1}(p)$. $\hat{\theta}_1 = [\hat{\rho}_1, \hat{\eta}_1]^T$ where $C/D = 1$.

### Step 2 Iterative iv estimation do

(i) iter = iter + 1

Generate the instrumental vector $y^{iv}(t)$ from (39) and from the estimated parameters $\hat{\rho}_{\mathrm{iter}-1}$.

(ii) Update the CT filters $F_{c,\mathrm{iter}}(p)$ in (37) with the new estimated parameter vector $\hat{\rho}_{\mathrm{iter}-1}$. Then, evaluate the filtered derivatives of the input $u(t)$, output $y^*(t)$ and instruments $y^{iv}(t)$.

(iii) Estimate the parameter vector of the noise model $\hat{\eta}_{\mathrm{iter}}$ by using an ARMA estimation algorithm on the output error $\hat{\xi}(t_k)$ (43).

(iv) Update the DT filters $F_{d,\mathrm{iter}}(q^{-1})$ in (38) with the estimated parameter vector $\hat{\eta}_{\mathrm{iter}}$. Then, apply the DT filter to the derivatives of signals calculated in (ii), namely the derivatives of $u_{f_c}(t_k)$, $y_{f_c}^*(t_k)$ and $y_{f_c}^{iv}(t_k)$, in order to define the regression vector $\varphi_f^*(t_k)$ and the instrumental vector $\varphi_f^{iv}(t_k)$, as in (40).

(v) From the filtered signals, compute the new parameter vector $\hat{\rho}_{\mathrm{iter}}$ according to (41). The global parameter vector $\hat{\theta}_{\mathrm{iter}}^{rivcf} = [\hat{\rho}_{\mathrm{iter}}, \hat{\eta}_{\mathrm{iter}}]^T$ is then obtained at iteration iter.

while $\max_j |\hat{\theta}_{\mathrm{iter}}^{j,rivcf} - \hat{\theta}_{\mathrm{iter}-1}^{j,rivcf}| > \epsilon$

where $\hat{\theta}_{\mathrm{iter}}^{j,rivcf}$ corresponds to the $j^{th}$ element of the parameter vector $\hat{\theta}_{\mathrm{iter}}^{rivcf}$ obtained at iteration iter.

### Step 3 Parametric error estimation

Compute the covariance matrix of the parametric error of the estimates $\hat{\rho}^{rivcf}$ with (44).

## 3. NUMERICAL EXAMPLE

A simulation example illustrates the performance of the *rivcf* algorithm and of the proposed full *rivcf* method, in case of colored output noise, in system identification using fractional models. The simulated model is given by:

$$ G_0(p) = \frac{K(-Tp^{\nu} + 1)}{((\frac{p}{\omega_1})^{2\nu} + 2\zeta_1(\frac{p}{\omega_1})^{\nu} + 1) ((\frac{p}{\omega_2})^{2\nu} + 2\zeta_2(\frac{p}{\omega_2})^{\nu} + 1)}, \quad (45) $$

with $\nu = 0.5$, $K = -1$, $T = 0.5$, $\omega_1 = 0.2 \text{ rad/s}$, $\zeta_1 = -0.4$, $\omega_2 = 1 \text{ rad/s}$ and $\zeta_2 = -0.65$.

This fractional system has two resonant modes: one at $\omega_1 = 0.2 \text{ rad/s}$ with $\zeta_1 = -0.4$ and another one at $\omega_2 = 1 \text{ rad/s}$ with $\zeta_2 = -0.65$ (for stable fractional systems, a resonant peak may occur for negative pseudo-damping factors $\zeta_{1,2}^3$). Furthermore, the system has a positive $s^{\nu}$-zero at $s^{\nu} = 2$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Fig. a): Power spectral density of u
Fig. b) Noisy output $y^*$ and noise-free output y</figcaption></figure>

The input signal $u(t)$, a pseudo random binary sequence of maximum length is used to excite the system (45) over its whole bandwidth. Its amplitude switches between -5 and 5 and the power spectral density is given on Fig. 2-a). The sampling interval is fixed to $T_s = 5.10^{-2} s$. Based on the CT process (45), the following hybrid system is considered (see Fig. 2-b))

$$ \left\{ \begin{array}{l} y(t) = G_0(p) u(t) \\ y^*(t_k) = y(t_k) + \xi(t_k), \end{array} \right. \quad (46) $$

where $\xi(t_k)$ is a colored additive noise derived from an ARMA(2,1) process $1/C(q^{-1})$:

$$ \frac{C(q^{-1})}{D(q^{-1})} = \frac{1 + 0.92q^{-1}}{1 - 1.96q^{-1} + 0.97q^{-2}}, \quad (47) $$

with a signal-to-noise ratio of 10dB.

The instrumental variables are computed according to the iterative algorithm described in §2.4 in order to ripen the parametric estimation. The system model structure is assumed to be known and only the model parameters are estimated. The model is set to:

$$ G(p) = \frac{b_1 p^{\nu} + b_0}{a_4 p^{4\nu} + a_3 p^{3\nu} + a_1 p^{2\nu} + a_0 p^{\nu} + 1}. \quad (48) $$

The commensurate order is set to the true value 0.5. However, the commensurate order is not always known; when unknown, the optimal commensurate order can be estimated. Moreover, the noise model structure is supposed known:

$$ \frac{C(q^{-1})}{D(q^{-1})} = \frac{1 + c_1 q^{-1}}{1 + d_1 q^{-1} + d_2 q^{-2}}. \quad (49) $$

A Monte-Carlo Simulation (MCS) of 200 runs is carried out, with a white noise input to the ARMA model being selected randomly at each run. The MCS results obtained from the *rivcf* and the full *rivcf* algorithms are shown in Table 1. From the *rivcf* columns of Table 1, it is shown that this algorithm delivers satisfactory results, however it is less accurate than the full *rivcf* method. Moreover, the *rivcf* algorithm clearly improves the statistical efficiency: in Table 1, $\sigma_{\theta_j}$ of *rivcf* method are much lower than the ones in *rivcf*. Furthermore, an iterative algorithm of the *rivcf* method clearly improves the statistical estimation.

<footnote n="3">3 When $\nu \neq 1$, $\zeta$ cannot be considered as a damping factor.</footnote>
<table><caption>Table 1. Monte-Carlo simulation results: parameters of the models (48) and (49) evaluated with the *srivcf*, *rivcf* (SR: Single Run) and *rivcf* methods ($\bar{\theta}$ is the mean and $\sigma$ is the standard deviation), with SNR = 10 dB</caption><thead><tr><th rowspan="2"></th><th colspan="3">srivcf</th><th colspan="3">rivcf (SR)</th></tr><tr><th>θ</th><th>θ̃</th><th>σ</th><th>θ̃</th><th>σ</th><th>θ̃</th><th>σ</th></tr></thead><tbody><tr><td>a<sub>4</sub></td><td>25</td><td>24.837</td><td>1.444</td><td>25.001</td><td>0.559</td><td>25.000</td><td>0.396</td></tr><tr><td>a<sub>3</sub></td><td>-36.5</td><td>-36.265</td><td>1.688</td><td>-36.487</td><td>0.732</td><td>-36.505</td><td>0.512</td></tr><tr><td>a<sub>2</sub></td><td>31.2</td><td>31.073</td><td>1.117</td><td>31.190</td><td>0.559</td><td>31.201</td><td>0.356</td></tr><tr><td>a<sub>1</sub></td><td>-5.3</td><td>-5.277</td><td>0.277</td><td>-5.301</td><td>0.097</td><td>-5.301</td><td>0.085</td></tr><tr><td>b<sub>1</sub></td><td>0.5</td><td>0.516</td><td>0.138</td><td>0.499</td><td>0.012</td><td>0.500</td><td>0.011</td></tr><tr><td>b<sub>0</sub></td><td>-1</td><td>-1.003</td><td>0.057</td><td>-0.998</td><td>0.023</td><td>-0.999</td><td>0.021</td></tr><tr><td>c<sub>1</sub></td><td>0.920</td><td></td><td></td><td>0.922</td><td>4.94 10<sup>-3</sup></td><td>0.920</td><td>4.70 10<sup>-3</sup></td></tr><tr><td>d<sub>2</sub></td><td>0.970</td><td></td><td></td><td>0.968</td><td>3.34 10<sup>-3</sup></td><td>0.970</td><td>3.06 10<sup>-3</sup></td></tr><tr><td>d<sub>1</sub></td><td>-1.960</td><td></td><td></td><td>-1.964</td><td>3.29 10<sup>-3</sup></td><td>-1.960</td><td>3.04 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. b)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Fig. a): Bode diagrams of the CT system (–), and its estimated models for 200 runs of Monte-Carlo with the *srivcf* (–.–) and the *rivcf* (– –) methods
Fig. b): Bode diagrams of the DT system (–), and its estimated models for 200 runs of Monte-Carlo with the *rivcf* (– –) methods</figcaption></figure>

In Fig. 3 -a), the 200 full *rivcf* identified models coincide with the true system; moreover, the 200 *srivcf* identified models also converge to the true system with a higher parameter variance. The true noise model and its 200 full *rivcf* identified noise models are plotted in the Bode diagram of Fig. 3 -b). As expected, the true system and the identified models coincide perfectly.

## 4. CONCLUSION

After a short introduction on fractional models and a synopsis of the *srivcf* method, the full *rivcf* algorithm for identifying hybrid Box-Jenkins fractional systems is presented when being in presence of colored noise. Even though the *srivcf* method produces satisfactory results, it is shown through Monte-Carlo simulations that the full *rivcf* method improves clearly the statistical efficiency.

When differentiation orders are unknown, as in practical physical systems, they can be estimated using a gradient-

based approach (see Victor et al. (2009) for the white noise case).

## REFERENCES

O. Cois, A. Oustaloup, T. Poinot, and J.-L. Battaglia. Fractional state variable filter for system identification by fractional model. In *6th European Control Conference ECC'01*, Porto, Portugal, 2001.

H. Garnier and L. Wang. *Identification of CT models from sampled data*. Springer-Verlag, 2008.

L. Ljung. *System identification – Theory for the user*. Prentice-Hall, Up. Sad. River, N.J., USA, 2 ed., 1999.

R. Malti, S. Victor, A. Oustaloup, and H. Garnier. An optimal IV method for CT fractional model identification. In *The 17th IFAC World Congress (IFAC'08)*, Seoul, Korea, July 2008.

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

K.B. Oldham and J. Spanier. *The fractionalal calculus – Theory and Applications of Differentiation and Integration to Arbitrary Order*. Academic Press, N.-Y., 1974.

A. Oustaloup. *La dérivation non-entière*. Hermès - Paris, 1995.

S.G. Samko, A.A. Kilbas, and O.I. Marichev. *Fractional integrals and derivatives: theory and applications*. Gordon and Breach Science, 1993.

T. Söderström and P. Stoica. *System Identification. Series in Systems and Control Engineering*. P. H., 1989.

S. Victor, R. Malti, and A. Oustaloup. IV method with optimal fractional differentiation order for CT system identification. In *15th IFAC Symposium on System Identification (SYSID'2009)*, Saint Malo, France, July 2009.

P.C. Young. Optimal IV identification and estimation of CT TF models. In *15th World IFAC Congress, Barcelona (Spain)*, 2002. IFAC, Elsevier.

P.C. Young and A.J. Jakeman. Refined IV methods of time-series analysis: Part III, extensions. *International Journal of Control*, 31:741–764, 1980.

P.C. Young, H. Garnier, and M. Gilson. An optimal IV approach for identifying hybrid CT Box-Jenkins models. In *14th IFAC Symposium on System Identification (SYSID'2006)*, pages 225–230, Newcastle, Australia, 2006.