# Model order identification for fractional models

Stéphane Victor, Rachid Malti

► To cite this version:

Stéphane Victor, Rachid Malti. Model order identification for fractional models. European Control Conference ECC'13, Jul 2013, Zurich, Switzerland. 10.23919/ECC.2013.6669589 . hal-00804797

HAL Id: hal-00804797

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

Submitted on 29 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
# Model order identification for fractional models

Stéphane Victor, Rachid Malti

Université de Bordeaux, IMS, CNRS UMR 5218,
351 cours de la Libération, 33405 Talence Cedex, France
{stephane.victor, rachid.malti}@ims-bordeaux.fr

## Abstract

This paper deals with continuous-time system identification using fractional differentiation models. So far, no algorithm exists concerning model order identification of fractional models. The “simplified” refined instrumental variable method is proposed to estimate parameters of fractional differential equation models when all the fractional orders are assumed known. Then an optimization approach based on this instrumental variable estimator is presented. This two-stage algorithm, called *oosrivcf* is then used for model order identification by using simultaneously two criteria: the Young Information Criterion (YIC) and the $R_T^2$ criterion.

## I. INTRODUCTION

Fractional calculus dates back to Liouville and Riemann in the nineteenth century and excels in modeling diffusive phenomena. Oldham and Spanier [1] expressed diffusive phenomena in semi-infinite media where diffusive systems are modeled with transfer functions involving differentiation orders multiples of 0.5. In electrochemistry, it is proven that the diffusion of charges in acid batteries is governed by Randles models [2], [3] that involve a half order integrator. In semi-infinite thermal systems, Battaglia et al. [4] have shown that the exact solution of heat equation links thermal flux to a half order derivative of surface temperature on which the flux is applied. Foucault currents inside rotor bars in induction machines obey to diffusive phenomena modeled by fractional models [5].

Instrumental variable methods for system identification have been developed since the seventies. For a detailed description refer to [6], [7], [8] and all references therein. **Refined Instrumental Variables for Continuous system (rivc)**, and its simplified version (*srivc*) when the additive measurement noise is white, were developed in [6], [9]. When the model lies in the same class as the true system, the *rivc* method and consequently the *srivc* can be interpreted in optimal statistical terms yielding consistent estimates with minimum variance.

In system identification with rational models, where only the coefficients are estimated, the model order remains unchanged, whereas estimating coefficients and differentiation orders in fractional models is not a trivial task given that estimating the differentiation orders in an iterative algorithm changes the model order at each iteration. In this paper, the *srivc* algorithm helps for linear coefficient estimation, and a gradient-based approach allows the differentiation order estimation. [10] have exposed this two-stage algorithm with an application to a thermal system. However, to the best of authors' knowledge, no algorithm has ever been proposed in the literature to fix the model structure, and consequently to fix optimally the number of parameters.

Fractional models allow modeling infinite dimensional systems with reduced number of parameters. Since it is always possible to find a high order rational model equivalent to a low dimension fractional model [11], the use of high dimension fractional models is of limited interest.

## Mathematical background

Differentiation to an arbitrary order $\gamma$, with $\gamma \in \mathbb{R}$ was defined by Riemann and Liouville as being an integer derivative of order $\lfloor \gamma \rfloor + 1$ ($\lfloor \cdot \rfloor$ stands for the floor operator) of a non-integer integral $\mathbf{I}$ [12]:

$$p^{\gamma}x(t) = p^{\lfloor \gamma \rfloor + 1} \left( \mathbf{I}^{\lfloor \gamma \rfloor + 1 - \gamma} x(t) \right) \triangleq \left( \frac{d}{dt} \right)^{\lfloor \gamma \rfloor + 1} \left( \frac{1}{\Gamma(\lfloor \gamma \rfloor + 1 - \gamma)} \int_{0}^{t} \frac{x(\tau) d\tau}{(t-\tau)^{\gamma-\lfloor \gamma \rfloor}} \right), \quad (1)$$

where $t > 0$, $\gamma \in \mathbb{R}_+^*$. The Euler's $\Gamma$ function is defined in the set of strictly positive real numbers $x \in \mathbb{R}_+^*$ as:

$$\Gamma(x) = \int_{0}^{\infty} e^{-t} t^{x-1} dt, \qquad (2)$$

and extended to the rest of the real numbers excepting negative integers $x \in \mathbb{R}^* \setminus \mathbb{N}^-$ by analytical continuity.
It is assumed that the system input $u(t)$ and the noise-free system output $y(t)$ are related by the following fractional differential equation,

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

where $p$ stands for the differentiation operator $p = \frac{d}{dt}$, and so $p^{\mu}u(t)$ denotes the $\mu$-th time-derivative of the continuous-time signal $u(t)$, $(a_j, b_i) \in \mathbb{R}^2$, and the differentiation orders

$$0 < \alpha_1 < \dots < \alpha_N, \quad 0 \le \beta_0 < \beta_1 < \dots < \beta_M, \qquad (4)$$

are ordered for identifiability purposes.

It is easy to define a symbolic representation of the fractional dynamic system governed by (3) using the transfer function:

$$G(p) = \frac{B(p)}{A(p)} = \frac{\sum_{i=0}^{M} b_i p^{\beta_i}}{1 + \sum_{j=1}^{N} a_j p^{\alpha_j}}. \qquad (5)$$

$G(p)$ is said commensurate of order $\gamma$, when rewritten as:

$$G(p) = \frac{\sum_{i=0}^{m} \tilde{b}_i p^{i\gamma}}{1 + \sum_{j=1}^{n} \tilde{a}_j p^{j\gamma}}, \qquad (6)$$

where $m = \frac{\beta_M}{\gamma}$ and $n = \frac{\alpha_N}{\gamma}$ are integers and $\forall i' \in \{0, 1, \dots, m\}, \forall j' \in \{1, \dots, n\}$:

$$\left\{ \begin{array}{ll} \tilde{b}_{i'} = b_i & \text{if } \exists i \in \{0, 1, \dots, m\} \text{ such that } i'\gamma = \beta_i \\ \tilde{b}_{i'} = 0 & \text{otherwise} \\ \tilde{a}_{j'} = a_j & \text{if } \exists j \in \{1, \dots, n\} \text{ such that } j'\gamma = \alpha_j \\ \tilde{a}_{j'} = 0 & \text{otherwise.} \end{array} \right. \qquad (7)$$

In rational transfer functions, $\gamma$ equals 1 and usually numerator $\alpha_N$ and denominator $\beta_M$ orders are both fixed, then all coefficients $b_i, i = 0, \dots, M$ and $a_j, j = 1, \dots, N$ are estimated. Generally, no care is taken to check whether any intermediate coefficient, as in (7), equals zero.

Time-domain simulation of fractional systems is an extensively studied topic in the literature [11], [13], [14]. System identification algorithms proposed in this paper may be used with any time-domain simulation algorithm of fractional systems. It is however definitely of his responsibility to correctly simulate the fractional system with negligible simulation errors, in order to be able to consistently estimate the fractional model parameters. Otherwise, simulation errors would affect the estimated parameters. The simulation algorithm used in this paper is based on [11] approximation of fractional operator in a desired frequency band by a recursive distribution of zeros and poles.

**Stability theorem ([15], extended):** A commensurate transfer function, with a commensurate order $\gamma$ as in (6), $G(p) = \frac{T(p^\gamma)}{R(p^\gamma)}$, where $T$ and $R$ are coprime polynomials, is BIBO stable if and only if:

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

and

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

The stability conditions are used later when estimating the commensurate order, which must be in the interval (8) for stable systems.

## II. PROBLEM FORMULATION

For simplicity, the formulation and solution of the estimation problem is restricted to the case of linear, single-input, single-output system, the analysis being extendable to multi-variable fractional systems (see e.g.[16]).

It is assumed that the input $u(t)$ and the noise-free output $y(t)$ are related by the constant coefficient differential equation (3), which can also be written in the compact model form (5), where $B(p)$ and $A(p)$ are assumed to be coprime, and the system asymptotically stable.

The data are collected at regular time-instants from $t = 0$ to $t = T_{final}$ (K samples), with a sampling period $h$ small enough, so that the approximation errors while computing numerically the fractional derivatives are negligible. The number
of samples is assumed to be large enough to guarantee convergence of the estimated parameters to the true ones. The quasi-stationary input signal {$u(t), 0 \le t \le T_{final}$} applied to the system is persistently exciting, and gives rise to an output signal {$y(t), 0 \le t \le T_{final}$}. The noise-free output y is supposed to be corrupted by an additive white measurement noise $\xi$, normally distributed with a zero mean and $\sigma^2$ variance, considered at discrete instants. The complete equation can be written in the form:

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

where $y^*(t_k)$ is the sampled value of the unobserved noisy output $y(t)$. Given the discrete-time sampled nature of the data, the usual assumption is that a discrete-time noise $\xi(t_k)$ is associated with the sampled data (see e.g. [17]). It is easier to consider white noise in the discrete-time case, since the concept of continuous-time white noise is a non-trivial extension of the discrete-time case.

When the model (5) is used, m Estimating all differentiation orders requires to determine $N + M + 1$ terms along with the $N + M + 1$ transfer function coefficients. Moreover, if $N$ and/or $M$ are high, the number of local minima may increase and the nonlinear optimization algorithms, with high complexity, may fail to converge to the global minimum. To reduce the number of parameters, a commensurate model (6) can be chosen, where only one differentiation order and hence $N + M + 2$ parameters are estimated. The parameter vector $\theta$ in (??) is then reduced to the $N + M + 1$ transfer function coefficients as in (??) added to a single commensurate differentiation order:

$$ \mu = \gamma. \qquad (11) $$

In this case, the searching domain of the commensurate order is held for $\gamma \in (0, 2)$, which corresponds to the stability domain as specified in the extended Matignon's theorem.

### III. MODEL ORDER IDENTIFICATION FOR CONTINUOUS-TIME FRACTIONAL MODELS

#### A. Simplified refined instrumental variable for continuous-time fractional models

In this section all differentiation orders are assumed to be known and fractional transfer function coefficients are estimated. The *srivc* [9], [18] algorithm derives from the *sriv* [19], [6] algorithm for discrete-time transfer function model identification and have been successfully used for more than three decades in various contexts. This was evolved by converting the maximum likelihood (ML) estimation equations to a pseudo-linear form involving optimal prefilters (see e.g. [7]). The iterative *srivcf* algorithm (*srivc* for fractional models), known to deliver optimal estimates when the additive measurement noise is white, is recalled (see [20], [10] for the details). Moreover, it has recently been established in the rational continuous-time case in [21, theorem 1] that the *srivc* algorithm converges asymptotically in one iteration under some mild assumptions<fnref n="1" />. In the fractional case, the convergence of the *srivcf* algorithm follows exactly the same scheme under the same assumptions.

All linear transformations can be interpreted as lowpass filtering. Not only direct differentiation of noised signals is avoided, but also it allows to specify the frequency band into which the model and system adequacy is sought. Moreover the filter frequency response should match its system transfer function. The *srivcf* uses a filter updated iteratively which avoids choosing a *priori* its design parameters.

As in the rational case, if the quasi-stationary input signal is persistently exciting and the noise signal is white with zero-mean, the *srivcf* estimator is asymptotically unbiased thanks to the *iv* mechanism. When the system is in the true model class, the *iv* estimate provides a consistent estimate under the following two conditions

$$ \left\{ \begin{array}{l} E \left[ \varphi_{iv,f}^{\text{iter}}(t_k) \varphi_f^{\text{T}}(t_k) \right] \text{ is non-singular,} \\ E \left[ \varphi_{iv,f}^{\text{iter}}(t_k) \xi(t_k) \right] = 0. \end{array} \right. \qquad (12) $$

For lack of space, only the main steps of the *srivcf* algorithm are recalled as follows. The interested reader may refer to [10] for more details.

**Step 1 Initialization:**

iter=0

Generate an initial<fnref n="2" /> estimate of coefficient parameter vector $\hat{\rho}^0$.

**Step 2 Iterative estimation:**

**do**

(i) Generate the auxiliary model output

$$ y_{iv}^{\text{iter}}(t) = \frac{\hat{B}(p, \hat{\rho}^{\text{iter}})}{\hat{A}(p, \hat{\rho}^{\text{iter}})} u(t). \qquad (13) $$

<footnote n="1">¹as qualified by the authors of the paper.</footnote>

<footnote n="2">²Least-squares or instrumental variable methods combined with fractional state variable filters can be used (see [20]).</footnote>
(ii) Update the filter

$$F_{\sigma}^{\text{iter}}(p) = \frac{p^{\sigma}}{\hat{A}(p, \hat{\rho}^{\text{iter}})}, \qquad (14)$$

and compute the filtered derivatives of the input, the output, and the instruments:

$$\begin{aligned}
p^{\beta_i} u_f(t) &= F_{\beta_i}^{\text{iter}}(p)u(t), \\
p^{\alpha_j} y_f^*(t) &= F_{\alpha_j}^{\text{iter}}(p)y^*(t), \\
p^{\alpha_j} y_{\text{iv},f}^{\text{iter}}(t) &= F_{\alpha_j}^{\text{iter}}(p)y_{\text{iv}}^{\text{iter}}(t).
\end{aligned}$$

(iii) Generate the regression $\varphi_f(t)$ and the instrumental variable $\varphi_{\text{iv},f}^{\text{iter}}(t)$ vectors :

$$\begin{aligned}
\varphi_f(t) &= \begin{bmatrix} p^{\beta_0} u_f(t), \dots, p^{\beta_M} u_f(t), \\ -p^{\alpha_1} y_f^*(t), \dots, -p^{\alpha_N} y_f^*(t) \end{bmatrix}^T \\
\varphi_{\text{iv},f}^{\text{iter}}(t) &= \begin{bmatrix} p^{\beta_0} u_f(t), \dots, p^{\beta_M} u_f(t), \\ -p^{\alpha_1} y_{\text{iv},f}^{\text{iter}}(t), \dots, -p^{\alpha_N} y_{\text{iv},f}^{\text{iter}}(t) \end{bmatrix}^T
\end{aligned}$$

(iv) Based on these prefiltered data, compute the new estimates $\hat{\theta}_{\text{iter}}^{\text{srivcf}}$

$$\hat{\rho}^{\text{iter}+1} = (\mathbf{\Phi}_{\text{iv},f}^{\text{iter}} \mathbf{\Phi}_f^{\text{T}})^{-1} \mathbf{\Phi}_{\text{iv},f}^{\text{iter}} \mathbf{Y}_f^*, \qquad (15)$$

with

$$\begin{aligned}
\mathbf{\Phi}_{\text{iv},f}^{\text{iter}} &= [\varphi_{\text{iv},f}^{\text{iter}}(0) \quad \dots \quad \varphi_{\text{iv},f}^{\text{iter}}(T_{\text{final}})]^{\text{T}}, \\
\mathbf{\Phi}_f &= [\varphi_f(0) \quad \dots \quad \varphi_f(T_{\text{final}})]^{\text{T}}, \\
\mathbf{Y}_f^* &= [y_f^*(0) \quad \dots \quad y_f^*(T_{\text{final}})]^{\text{T}}.
\end{aligned}$$

$$\textbf{while } \max_j \left| \frac{\hat{\rho}_j^{\text{iter}+1} - \hat{\rho}_j^{\text{iter}}}{\rho_j^{\text{iter}+1}} \right| > \epsilon_1,$$

where $\hat{\rho}_j^{\text{iter}}$ is the j-th element of the estimated parameter vector at the iteration iter.

Step 3 **Parametric error estimation**

Compute the estimated parametric error covariance matrix

$$\hat{\mathbf{P}}_{\rho} = \hat{\sigma}^2 (\mathbf{\Phi}_{\text{iv},f} \mathbf{\Phi}_{\text{iv},f}^{\text{T}})^{-1}, \qquad (16)$$

where $\hat{\sigma}^2$ is the empirical estimation of the noise variance, and $\mathbf{\Phi}_{\text{iv},f}$ corresponds to $\mathbf{\Phi}_{\text{iv},f}^{\text{iter}}$ computed at the last iteration.

## B. Differentiation order estimation

When differentiation orders are unknown, as it is often the case in practice, it is helpful to consider order estimation along with transfer function coefficient estimation. An algorithm is proposed to estimate the parameters in two stages. This algorithm, named order-optimization-srivcf or oosrivcf for short, uses the srivcf algorithm for coefficient estimation and is combined with a gradient-based algorithm for differentiation order optimization. Two variants of order optimization algorithm are recalled [10]. Either all differentiation orders are set as integral multiples of a commensurate order, as in the model (??), and the best commensurate order is computed in the interval (0, 2), or all differentiation orders, as in (5), are estimated. The former variant is used as a good initial hit for the latter variant. The estimation problem is formulated as a minimization problem of the $\ell_2$-norm:

$$J(\theta) = \frac{1}{2} \|\varepsilon(t, \theta)\|^2, \qquad (17)$$

of the output error

$$\varepsilon(t, \theta) = y^*(t) - G(p, \theta)u(t), \qquad (18)$$

with respect to the vector $\mu$, as in (??) or (11).

A Gauss-Newton algorithm (see e.g. [22]) is used for the estimation of all differentiation orders $\mu$ as defined in (??). For lack of space, only the main steps of the oosrivcf algorithm are recalled as follows:

Step 1 **Initialization**

iter = 0
Initialize $\mu^0$ and compute $\rho^0$ with the srivcf method.
From $\theta^0 = [\rho^0, \mu^0]$, compute $J(\theta^0)$.

## Step 2 Gauss-Newton optimization
do

Initialize $\lambda = \Lambda$ (usually to 1)
do

(i) Compute the error sensitivity function:

$$ \frac{\partial \varepsilon}{\partial \mu} = \left[ -\frac{\partial \hat{y}}{\partial \beta_0}, \dots, -\frac{\partial \hat{y}}{\partial \beta_M}, -\frac{\partial \hat{y}}{\partial \alpha_1}, \dots, -\frac{\partial \hat{y}}{\partial \alpha_N} \right]^T, \quad (19) $$

where the output sensitivity functions are computed numerically.

(ii) Compute the gradient

$$ \frac{\partial J}{\partial \mu} = \frac{\partial \varepsilon^T}{\partial \mu} \varepsilon, \quad (20) $$

and the approximated Hessian

$$ H = \frac{\partial \varepsilon^T}{\partial \mu} \frac{\partial \varepsilon}{\partial \mu}. \quad (21) $$

(iii) Compute the new order estimate:

$$ \mu^{\text{iter}+1} = \mu^{\text{iter}} - \lambda \left[ H^{-1} \frac{\partial J}{\partial \mu} \right]_{\mu^{\text{iter}}}. $$

(iv) Compute $\rho^{\text{iter}}$ using the srivcf.

(v) Evaluate the error criterion: $J(\theta^{\text{iter}+1})$.

(vi) $\lambda = \lambda/2$

while $J(\theta^{\text{iter}+1}) > J(\theta^{\text{iter}})$

iter = iter + 1

while $\max_l \left| \frac{\mu_l^{\text{iter}} - \mu_l^{\text{iter}-1}}{\mu_l^{\text{iter}}} \right| > \epsilon_2$

where $\mu_l^{\text{iter}}$ corresponds to $l$-th element of the order vector $\mu^{\text{iter}}$ at iteration iter. $l$ equals one in case of commensurate order estimation.

## Step 3 Parametric error estimation

After convergence, the covariance matrix $\hat{\mathbf{P}}_{\theta}$ associated to the estimate $\hat{\theta}$ can be estimated [22, Theorem 9.1] by:

$$ \hat{\mathbf{P}}_{\theta} = \hat{\sigma}^2 \mathcal{H}^{-1}, \quad (22) $$

where $\hat{\sigma}^2$ is, as previously, the empirical estimate of noise variance and $\mathcal{H}$ is the approximate Hessian computed towards all the estimated parameters:

$$ \mathcal{H} = \frac{\partial \varepsilon^T}{\partial \theta} \frac{\partial \varepsilon}{\partial \theta} = \frac{\partial \varepsilon^T}{\partial [\rho^T \mu^T]^T} \frac{\partial \varepsilon}{\partial [\rho^T \mu^T]^T}, \quad (23) $$

with $\frac{\partial \varepsilon}{\partial \mu}$ defined in (19) and

$$ \frac{\partial \varepsilon}{\partial \rho} = - \left[ \frac{\partial \hat{y}}{\partial b_0}, \dots, \frac{\partial \hat{y}}{\partial b_M}, \frac{\partial \hat{y}}{\partial a_1}, \dots, \frac{\partial \hat{y}}{\partial a_N} \right]^T, \quad (24) $$

where $\forall r = 0, \dots, M$, and $\forall k = 1, \dots, N$

$$ \frac{\partial \hat{y}}{\partial b_r} = \frac{\hat{p}^{\hat{\beta}_r}}{1 + \sum_{j=1}^{N} \hat{a}_j p^{\hat{\alpha}_j}} u(t), \quad (25) $$

$$ \frac{\partial \hat{y}}{\partial a_k} = - \frac{\sum_{i=0}^{M} \hat{b}_i p^{\hat{\beta}_i + \hat{\alpha}_k}}{\left( 1 + \sum_{j=1}^{N} \hat{a}_j p^{\hat{\alpha}_j} \right)^2} u(t). \quad (26) $$
Remarks:

1) In the Gauss-Newton method, the advantage of computing the approximated Hessian is that the second derivative of ε is not required; hence it is less time consuming. Moreover, it allows to obtain always a positive semi-definite matrix which is a necessary convergence condition. If the criterion is convex, the algorithm reaches the global minimum. Otherwise, only the convergence to a local minimum is guaranteed. The λ parameter is used to adjust the step in case of slow convergence or oscillations around a minimum.

2) The output sensitivity functions ($\forall r = 0, \dots, M$ and $\forall k = 1, \dots, N$)

$$ \frac{\partial \hat{y}}{\partial \beta_r} = \ln(p) \frac{\hat{b}_r p^{\hat{\beta}_r}}{1 + \sum_{j=1}^{N} \hat{a}_j p^{\hat{\alpha}_j}} u(t), \qquad (27) $$

$$ \frac{\partial \hat{y}}{\partial \alpha_k} = \ln(p) \hat{a}_k p^{\hat{\alpha}_k} \frac{\sum_{i=0}^{M} \hat{b}_i p^{\hat{\beta}_i}}{\left(1 + \sum_{j=1}^{N} \hat{a}_j p^{\hat{\alpha}_j}\right)^2} u(t). \qquad (28) $$

depend on $\ln(p)$ (log of the differential operator) which is not trivial to simulate in the time-domain. That is why they are computed numerically, using for example the central difference method.

3) The oosrivcf algorithm works similarly when considering the commensurate order; $\mu$ is then defined as in (11). Instead of using (19), the following error sensitivity function will be used:

$$ \frac{\partial \varepsilon}{\partial \mu} = \frac{\partial \varepsilon}{\partial \gamma} = -\frac{\partial \hat{y}}{\partial \gamma}, \qquad (29) $$

where

$$ \frac{\partial \hat{y}}{\partial \gamma} = \ln(p) \frac{\sum_{i=0}^{M} i \hat{b}_i p^{i \hat{\gamma}} + \sum_{i=0}^{M} \sum_{j=1}^{N} (i-j) \hat{b}_i \hat{a}_j p^{(i+j) \hat{\gamma}}}{\left(1 + \sum_{j=1}^{N} \hat{a}_j p^{j \hat{\gamma}}\right)^2} u(t) \qquad (30) $$

depends on $\ln(p)$. As previously, it will also be computed numerically.

In case all differentiation orders are estimated, it is highly recommended to estimate a commensurate model, and then to release the commensurability constraint in order to estimate all differentiation orders. The commensurate model constitutes generally a good initial hit when estimating all differentiation orders. Monte Carlo simulation results have been successfully carried out and presented in [10].

C. Model Order estimation

In the identification procedure, a key point involves the model order structure selection: the number of parameters $N$ and $M$ in the transfer function (5) need to be fixed. This point is named as Model order selection in rational systems. However, in fractional systems, model order $\alpha_N$ in (5) is estimated along with the coefficients. Hence, it is preferable to talk about structure selection. The method developed for rational SISO systems [9] is extended to fractional systems. Two statistical measures are computed and used to choose between a range of orders $M$ and $N$, which are defined as:

$$ R_T^2 = 1 - \frac{\hat{\sigma}_{\xi}^2}{\hat{\sigma}_y^2}, \qquad (31) $$

$$ YIC = \log \left\{ \frac{\hat{\sigma}_{\xi}^2}{\hat{\sigma}_y^2} \right\} + \log \left\{ \frac{1}{n_{\theta}} \sum_{i=1}^{n_{\theta}} \frac{\hat{\sigma}_{\xi}^2 \hat{p}_{ii}}{\hat{\theta}_i^2} \right\}, \qquad (32) $$

where $\hat{\sigma}_{\xi}^2$ and $\hat{\sigma}_y^2$ denote respectively the variance of the simulation error and the variance of the measured output, $\hat{p}_{ii}$ is the $i$-th diagonal element of the estimated parametric covariance matrix given in (22) and $\theta_i^2$ is the square of the $i$-th parameter estimate in the $\theta$ vector so that the ratio $\hat{p}_{ii}/\theta_i^2$ is a normalized measure of the uncertainty on the $i$-th parameter estimate.

$R_T^2$ is recognized as the coefficient of determination based on the simulation error. It is a normalized measure of how much of the output variance is explained by the deterministic system part of the estimated model and it will be close to 1 in low noise situations. However, this measure is not sufficient to avoid over-parametrization and identify a parsimonious model. The Young's information criterion YIC is more complex and provides a measure of how well parameters are defined statistically: the more negative YIC, the better the definition. However, it may lead to underestimate the orders $M$ and $N$. Both criteria are inspected to find the orders for which $R_T^2$ is sufficiently high to indicate a good explanation of the data and the YIC is sufficiently negative to indicate well-defined parameter estimates.
<table><caption>TABLE I
BEST MODEL STRUCTURES ESTIMATED WITH THE *oosrivcf* ALGORITHM ACCORDING TO YIC, AIC AND QUADRATIC CRITERIA.</caption><thead><tr><th rowspan="2">N</th><th rowspan="2">M</th><th colspan="2">YIC criterion</th><th>AIC criterion</th><th>quadratic criterion</th></tr><tr><th>J<sub>YIC</sub> (dB)</th><th>R<sup>2</sup><sub>T</sub></th><th>J<sub>AIC</sub> (dB)</th><th>J (dB)</th></tr></thead><tbody><tr><td>2</td><td>1</td><td>-6.886</td><td>0.802</td><td>-1.293</td><td>-7.031</td></tr><tr><td>3</td><td>1</td><td>-8.020</td><td>0.909</td><td>-1.632</td><td>-10.419</td></tr><tr><td>3</td><td>2</td><td>-4.021</td><td>0.909</td><td>-1.632</td><td>-10.420</td></tr><tr><td>4</td><td>1</td><td>-7.451</td><td>0.908</td><td>-1.627</td><td>-10.368</td></tr><tr><td>4</td><td>2</td><td>-6.174</td><td>0.909</td><td>-1.632</td><td>-10.416</td></tr><tr><td>4</td><td>3</td><td>-4.582</td><td>0.909</td><td>-1.629</td><td>-10.393</td></tr><tr><td>5</td><td>1</td><td>-2.661</td><td>0.908</td><td>-1.627</td><td>-10.369</td></tr><tr><td>5</td><td>2</td><td>-5.031</td><td>0.909</td><td>-1.629</td><td>-10.386</td></tr><tr><td>5</td><td>3</td><td>-4.574</td><td>0.909</td><td>-1.631</td><td>-10.412</td></tr><tr><td>5</td><td>4</td><td>-4.052</td><td>0.909</td><td>-1.631</td><td>-10.408</td></tr></tbody></table>

### D. Simulation example

Consider the case of estimating transfer function coefficients and all differentiation orders of the following data-generating system:

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

The input signal is a PRBS (see Fig.1).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Input, noise-free and noisy outputs, used as identification data, zoomed between $t = 0$ and 200s, identification data going up to 600s (SNR = 10dB and $h = 0.02s$)</figcaption></figure>

The output $y(t)$ is corrupted at the sampling instants (10) by an additive Gaussian white noise $\xi(t_k)$ with zero mean and a SNR of 10dB. It is highly recommended that a commensurate model, as in (6), be estimated first and used as an initial hit for estimating all orders.

The model structure is chosen as in (5) with $\beta_0 = 0$. As initial guess, the coefficients are fixed to 1, and the differentiation orders are multiple of 1.3.

The estimation synthesis is presented in Table I, where the first two columns define the model structure. Here the right model class is defined for $N = 3$ and $M = 1$ (marked in yellow in Table I). When the model is chosen in the right model class, the *oosrivcf* converged to the true parameters (up to to the noise effect):

$$G_1(p) = \frac{0.502}{0.502p^{2.796} + 1.503p^{1.197} + 1} \quad (34)$$

As a comparison, the AIC criterion:

$$AIC = \log \left\{ \left( 1 + \frac{2n_{\theta}}{K} \right) \| \varepsilon(t, \hat{\theta}) \|^2 \right\}, \quad (35)$$
and the quadratic criteria have been computed and displayed respectively on the last columns of Table I. Both criteria are not adapted for choosing the most suited fractional model class. Indeed, with the AIC criterion, multiple model structure seem to suit (the lowest criterion value are marked in green in Table I), and the quadratic criterion do not give the right model structure (marked in cyan in Table I).

The YIC and $R_T^2$ criteria have the advantage to directly give the most suited model structure (marked in red in Table I): the YIC should be the lowest and $R_T^2$ the highest.

## IV. CONCLUSION

The well-known *srive* algorithm extended to fractional models is used to estimate transfer function coefficients and it is combined to a gradient-based algorithm for fractional differentiation order estimation. This algorithm is presented as the *oosrивcf* algorithm, which combines simultaneously the transfer function coefficient and the differentiation order estimation. An initialization procedure is proposed consisting of estimating a commensurate model first, and then of using the obtained model as an initial hit for all order estimation. In order to get the best model structure, two criteria are required. The $R_T^2$ criterion, recognized as the coefficient of determination based on the simulate error is not sufficient to avoid over-parametrization and identify a parsimonious model. However, the Young's information criterion *YIC* provides a measure of how well parameters are defined statistically: the more negative *YIC*, the better the definition. It may lead to underestimate model orders. Both criteria are inspected to find the orders for which $R_T^2$ is sufficiently high to indicate a good explanation of the data and the *YIC* is sufficiently low to indicate well-defined parameter estimates. The performance of the proposed algorithms has been evaluated on a simulation example and compared to the AIC and quadratic criteria.

It will also be interesting to extend this study to deal with colored output noise by using hybrid Box-Jenkins models with continuous-time fractional input-output models and discrete-time noise models.

## REFERENCES

[1] K.B. Oldham and J. Spanier. *The fractionalal calculus - Theory and Applications of Differentiation and Integration to Arbitrary Order*. Academic Press, New-York and London, 1974.

[2] S. Rodrigues, N. Munichandraiah, and A.-K. Shukla. A review of state of charge indication of batteries by means of A.C. impedance measurements. *Journal of Power Sources*, 87:12–20, 2000.

[3] J. Sabatier, M. Aoun, A. Oustaloup, G. Grégoire, F. Ragot, and P. Roy. Fractional system identification for lead acid battery state charge estimation. *Signal processing*, 86(10):2645–2657, 2006.

[4] J.-L. Battaglia, O. Cois, L. Puigsegur, and A. Oustaloup. Solving an inverse heat conduction problem using a non-integer identified model. *Int. J. of Heat and Mass Transfer*, 44(14):2671–2680, 2001.

[5] A. Benchellal, S. Bachir, T. Poinot, and J.-C. Trigeassou. Identification of a non-integer model of induction machines. In *1st IFAC Workshop on Fractional Differentiation and its Applications (FDA)*, Bordeaux, France, 2004.

[6] P.C. Young. Parameter estimation for continuous-time models – a survey. *Automatica*, 17(1):23–29, 1981.

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

[8] T. Söderström and P. Stoica. *Instrumental variable methods for system identification*. Springer Verlag, New York, 1983.

[9] P.C. Young. Optimal IV identification and estimation of continuous-time TF models. In *15th World IFAC Congress, Barcelona (Spain)*, 2002. IFAC, Elsevier.

[10] S. Victor, R. Malti, H. Garnier, and A. Oustaloup. Parameter and differentiation order estimation in fractional models. *Automatica*, 49(4), 2013.

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

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

[13] M. Aoun, R. Malti, F. Lebron, and A. Oustaloup. Numerical simulations of fractional systems: an overview of existing methods and improvements. *An Int. J. of Nonlinear Dynamics and Chaos in Engineering Systems. Special issue: Fractional Derivatives and Their Applications.*, 38(1-4):117–131, 2004.

[14] B.T. Krishna. Studies on fractional order differentiators and integrators: A survey. *Signal Processing*, 91(3):386–426, 2011.

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

[16] A. J. Jakeman and P. C. Young. Refined instrumental variable methods of time-series analysis: Part ii, multivariable systems. *International Journal of Control*, 29:621–644, 1979.

[17] P.C. Young, H. Garnier, and M. Gilson. Identification of continuous-time models from sampled data. In H. Garnier and L. Wang, editors, *Identification of continuous-time models from sampled data*, chapter Refined Instrumental Variable Identification of CT Hybrid Box-Jenkins Models, pages 91–131. Springer-Verlag, 2008.

[18] H. Garnier and L. Wang. *Identification of continuous-time models from sampled data*. Springer-Verlag, 2008.

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

[20] R. Malti, S. Victor, and A. Oustaloup. Advances in system identification using fractional models. In ASME, editor, *Journal of Computational and Nonlinear Dynamics*, volume 3, pages 021401.1–021401.7, New York, USA, April 2008.

[21] X. Liu, J. Wang, and W.X. Zheng. Convergence analysis of refined instrumental variable method for continuous-time system identification. *IET Control Theory & Applications*, 5(7):868–877, May 2011.

[22] L. Ljung. *System identification – Theory for the user*. Prentice-Hall, Upper Saddle River, N.J., USA, 2 edition, 1999.