# Instrumental variable method with optimal fractional differentiation order for continuous-time system identification

Stéphane Victor, Rachid R. Malti, Alain Oustaloup

► To cite this version:

Stéphane Victor, Rachid R. Malti, Alain Oustaloup. Instrumental variable method with optimal fractional differentiation order for continuous-time system identification. 15th IFAC Symposium on System Identification (SYSID), Jul 2009, Saint-Malo, France. pp.904-909. hal-00401701

HAL Id: hal-00401701

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

Submitted on 8 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 method with optimal fractional differentiation order for continuous-time system identification

Stéphane Victor,\* Rachid Malti,\* Alain Oustaloup,\*

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

**Abstract:** this paper deals with continuous-time system identification using fractional differentiation models in a noisy output context. The simplified refined instrumental variable for continuous-time fractional systems (srivcf) is extended to optimize the commensurate order with a gradient-based method as the coosrivcf algorithm. Simulation analysis is used to demonstrate the performance of the proposed optimal instrumental variable scheme combined with order optimization.

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

## 1. INTRODUCTION

As detailed in the survey paper by Young (1981) and all references therein, Instrumental Variable (IV) approaches to system identification and parameter estimation have been developed since 1960s and 1970s. Recently, Refined Instrumental Variables for Continuous-system (Young and Jakeman (1980)) has been extended for identifying fractional continuous-time systems (Malti et al. (2008b)). The least-square (LS) method and the instrumental variable (IV) combined with state variable filter (SVF) are recalled on that paper. The model structure is fixed and only the coefficients were identified in a noisy context. To minimize the parameter variance, an iterative optimal instrumental variable method for fractional continuous-time model was implemented where the auxiliary model and the state variable filters are actualized at each iteration. The algorithm is initialized by a transfer function obtained from the least squares state variable filter estimates. Cois et al. (2001) are the first to extend the use of state variable filters to fractional systems. The authors suggest to use IV to eliminate bias. However, their instruments are not selected in an optimal way. The optimal approach involves a method of iterative-adaptive prefiltering based on a quasi-optimal statistical solution obtained from a given commensurate<fnref n="1" /> order. This optimal Instrumental Variables (IV) method presents the advantage of yielding asymptotically consistent estimates in the presence of measurement noise.

Thanks to the IV mechanism, it was shown in Malti et al. (2008b) that the srivcf method is asymptotically unbiased, provided the true system is in the right model class. To be in the same model class, a model needs to have the same dimension and the same commensurate differentiation order. The true dimension and the commensurate order are however usually unknown in practice which might induce

unmodeled dynamics. Consequently, it is proposed in this paper to fix the dimension of the system, as is usually done for rational models, and to combine the optimization of the commensurate order with a gradient-based method to the srivcf algorithm. The proposed method is independent from the way fractional differentiation and integration are simulated in the time-domain.

### 1.1 Mathematical background

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

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

where $(a_j, b_i) \in \mathbb{R}^2$, differentiation orders $\alpha_1 < \alpha_2 < \dots < \alpha_{m_A}$, $\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) $\gamma$, with $\gamma \in \mathbb{R}_+^*$ (set of strictly positive real numbers):

$$\mathbf{D}^{\gamma} \triangleq \left(\frac{d}{dt}\right)^{\gamma}$$

was defined in the 19th century by Riemann and Liouville. The $\gamma$-order fractional derivative of $x(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):

$$\mathbf{D}^{\gamma} x(t) = \mathbf{D}^{\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 (2)$$

where $t > 0$, $\gamma \in \mathbb{R}_+^*$, and the Euler's $\Gamma$ function is defined in the set of real numbers except zero and negative whole numbers ($x \in \mathbb{R}^* \setminus \mathbb{N}^-$) as:

<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>
$$ \Gamma(x) = \int_{0}^{\infty} e^{-t} t^{x-1} dt. \qquad (3) $$

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

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

This property allows to write the fractional differential equation (1), provided $u(t)$ and $y(t)$ equal 0 for all $t < 0$, in a transfer function form:

$$ F(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}}. \qquad (5) $$

Moreover if $F(s)$ is commensurate of order $\gamma$, then $F(s)$ can be rewritten as:

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

where $n = \frac{\beta_{m_B}}{\gamma}$ and $m = \frac{\alpha_{m_A}}{\gamma}$ are integers and:

$$ \left\{ \begin{array}{lll} \tilde{b}_i = b_i & \text{if } i\gamma = \beta_i & \text{and } \tilde{b}_i = 0 & \text{if } i\gamma \neq \beta_i \\ \tilde{a}_j = a_j & \text{if } j\gamma = \alpha_j & \text{and } \tilde{a}_j = 0 & \text{if } j\gamma \neq \alpha_j. \end{array} \right. \qquad (7) $$

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

Exactly the same principle will be applied here for fractional models: $\alpha_{m_A}$ and $\beta_{m_B}$ being fixed, all the $\tilde{b}_i$ and $\tilde{a}_j$ coefficients are estimated for a given commensurate order $\gamma_k$ which is estimated iteratively.

## 1.2 Time-domain simulation of fractional models

The main objective of this paper is to develop a method for time domain system identification able to work with any fractional model time domain simulation algorithm. There exists several simulation algorithms in the literature: the recursive approximation with a distribution of poles and zeros Oustaloup (1995), an integrator behavior in high and low frequencies Poinot and Trigeassou (2003) and more developments are detailed in Aoun et al. (2004). The objective being not to discuss the different simulation algorithms, hence, only the simulation algorithm used in the example of this paper is explained here.

Due to the consideration that real physical systems generally have bandlimited fractional behavior and due to the practical limitations of input and output signals (Nyquist frequency for the upper band and the spectrum of the input signal for the lower band), fractional operators are usually approximated by high order rational models. As a result, a fractional model and its rational approximation have the same dynamics within a limited frequency band. The most commonly used approximation of $s^\gamma$ in the frequency band $[\omega_A, \omega_B]$ is the recursive distribution of

<figure><img src="image_2.png" /><figcaption>Fig. 1. Recursive approximation of a fractional differentiator with poles and zeroes</figcaption></figure>

zeros and poles proposed by Oustaloup Oustaloup (1995) (see Fig. 1):

$$ s^{\nu} \rightarrow s_{[\omega_A, \omega_B]}^{\gamma} = C_0 \left( \frac{1 + \frac{s}{\omega_A}}{1 + \frac{s}{\omega_B}} \right)^{\gamma} \approx C_0 \prod_{k=1}^{N} \frac{1 + \frac{s}{\omega_k}}{1 + \frac{s}{\omega_k}} \qquad (8) $$

where $\omega_i = \alpha\omega'_i$, $\omega'_{i+1} = \eta\omega'_i$ and

$$ \gamma = 1 - \frac{\log \alpha}{\log \alpha\eta}, \qquad (9) $$

$\alpha$ and $\eta$ define the differentiation order $\gamma$. The bigger $N$ the better the approximation of the differentiator $s^\gamma$ in its frequency band. Moreover, it can be seen in Fig. 1 that the frequency bounds should also be well defined so that the approximation fits well the real fractional derivative $s^\gamma$ in a precise frequency band and so that the edge effects do not distort the approximation. The frequency band $[\omega_A, \omega_B]$ could be enlarged to get a better behavior.

## 2. SYSTEM IDENTIFICATION

To be in the right model class, an LTI model needs to have the same dimension and the same commensurate differentiation order as the system. The true dimension and the commensurate order are however usually unknown in practice, which might induce unmodeled dynamics. In the following, the dimension of the system is fixed, as is usually done in the rational case, and the optimization of the commensurate order is combined with the srivcf algorithm in a two stage procedure. First a commensurate order $\gamma_0$ is fixed. Then, for a fixed commensurate order, optimal coefficients are computed iteratively using the srivcf algorithm. Next, for the computed coefficients, an optimal commensurate order is calculated using nonlinear programming technics. This procedure is repeated until convergence as sketched in Fig. 2. Although such algorithms have no proof of convergence, they do often converge. The simplified refined IV (srivcf) is presented in the following section. The nonlinear programming algorithm for order optimization is detailed in 2.3.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Iterative order optimization algorithm: coosrivcf</figcaption></figure>

## 2.1 Simplified refined instrumental variable

The Refined Instrumental Variables for Continuous system (Young and Jakeman (1980)) (rivc) method and its simplified version srivc (when the additive measurement noise is assumed to be white) can be interpreted in optimal statistical terms (consistent estimators with minimum variance) when the system is in the right model class. Moreover, both methods provide an estimation of the covariance matrix and, therefore, an estimation on the parameter confidence bounds. Optimal instrumental variables are obtained from an auxiliary model which is iteratively updated. At the first iteration, the auxiliary model is initialized using a least squares estimate applied on state variable filters (SVF). The identification problem is first stated, then SVF method is briefly recalled and finally the srivcf algorithm is detailed.

Consider input $u(t)$ and output $y^*(t)$ data collected at regular samples $0, T_s, 2T_s, \dots, KT_s = T_{fin}$. The measured output is supposed to be corrupted by an additive measurement noise $\xi(t)$:

$$y^*(t) = y(t) + \xi(t), \quad (10)$$

where $y(t)$ is the hypothetical noise-free deterministic system output. The input $u(t)$ and the output $y(t)$ signals are supposed to be related by the fractional differential equation (1).

Equation error can be formulated as:

$$\varepsilon(t) = y^*(t) - \varphi^*(t)^T \theta, \quad (11)$$

where

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

and

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

## 2.2 Least square estimates

Minimizing the $L_2$ norm of $\varepsilon(t)$:

$$J = \int_{0}^{T_{fin}} (\varepsilon(t))^2 dt, \quad (14)$$

with respect to $\theta$, leads to the following least squares estimates:

$$\hat{\theta}_{LS} = \left[ \int_{0}^{T_{fin}} \varphi^*(t)^T \varphi^*(t) dt \right]^{-1} \int_{0}^{T_{fin}} \varphi^*(t)^T y^*(t) dt. \quad (15)$$

Or, after a numerical discretization, by defining $\mathbf{Y}$ as a column vector of the system output and $\mathbf{\Phi}$ as a regression matrix where the columns are fractional derivatives of the input and the output signals:

$$\mathbf{\Phi}^* = [\varphi^*(T_s) \ \varphi^*(2T_s) \ \dots \varphi^*(KT_s)]^T, \quad (16)$$

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

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

As in the rational case, direct fractional differentiations of noisy signals lead however to inaccurate results as it amplifies the noise. Cois et al. (2001) showed that the least square estimator is biased in presence of noisy output due to the correlation between $\varphi^*(t)$ and $y^*(t)$. Consequently, a linear transformation (low-pass filter) is applied to input and output signals: state variable filters (SVF) are introduced to the direct least squares estimates in order to filter the noise (see Cois et al. (2001) for the details or Malti et al. (2008b) with its combination with the instrumental variable which leads to unbiased estimates).

### 2.2.1. Srivcf estimates
As shown by Young (1981), the optimal estimator (asymptotically consistent with minimal variance) of the IV method is obtained when the SVF filters $F_{\gamma}^{opt}(s)$ are as:

$$F_{\gamma}^{opt}(s) = \frac{s^{\gamma}}{H_{o}(s)A(s)}, \quad (18)$$

where $H_o(s)$ is the real noise model, $A(s)$ the denominator of the transfer function of the real system; and the hypothetical noise free deterministic output $y(t)$ is used when evaluating the derivatives of the output. As a result, the optimal regression vector is obtained (by introducing the differential operator $p = \frac{d}{dt}$):

$$\varphi_f^{\text{opt}}(t)^T = F^{\text{opt}}(p) \varphi^*(t)^T = \begin{bmatrix} \mathbf{D}^{\beta_0} u_f(t) \cdots \mathbf{D}^{\beta_{m_B}} u_f(t) \\ -\mathbf{D}^{\alpha_1} y_f^*(t) \cdots -\mathbf{D}^{\alpha_{m_A}} y_f^*(t) \end{bmatrix}. \quad (19)$$

The optimal estimator requires to know many characteristics of the system to be identified. It will be obtained only if the real system and noise models are known; however the optimality is not achievable in practical cases.

A stochastic approach proposed by Young (1981, 2002) is the Simplified Refined Instrumental Variable for Continuous systems (srivc). This approach includes an adaptive pre-filtering based on the optimal solution when the noise $\xi$ is a white gaussian one ($H_o(s) = 1$). This method is a logical extension of the IV estimator with SVF. This estimator is extended to fractional models and is named srivcf (srivc for fractional models). A suitable error function is given by the output error:
$$e(t) = y(t) - \frac{B(p)}{A(p)}u(t) \quad (20)$$

The equation error can be rewritten:

$$\begin{aligned} e(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) \quad (21) \end{aligned}$$

By replacing $\frac{1}{A(p)}$ by $F(p)$ as the iteratively computed filter, the error becomes:

$$\begin{aligned} e(t) &= A(p)(F(p)y(t)) - B(p)(F(p)u(t)) \\ &= A(p)y_f(t) - B(p)u_f(t) \quad (22) \end{aligned}$$

To optimize the instruments, it is required to implement an iterative algorithm. Initialized by the $\hat{\theta}_{IVSVF}$ estimates for example (see Malti et al. (2008b)), the new derivatives of input and output signals are computed with the new optimal state variable filter:

$$F_{\gamma}^{(\text{iter})}(s, \hat{\theta}) = \frac{s^{\gamma}}{\hat{A}(s)} = \frac{s^{\gamma}}{1 + \sum_{j=1}^{m_A} \hat{a}_j s^{\alpha_j}}, \quad (23)$$

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

The instruments, the derivatives of the input and of the output, are computed at each iteration by:

$$\left\{ \begin{aligned} \mathbf{D}^{\beta_i} u_f(t) &= u(t) \otimes \mathcal{L}^{-1}\{F_{\beta_i}^{(\text{iter})}(s)\} \\ \mathbf{D}^{\alpha_j} y_f^*(t) &= y^*(t) \otimes \mathcal{L}^{-1}\{F_{\alpha_j}^{(\text{iter})}(s)\} \\ \mathbf{D}^{\alpha_j} y_f^{\text{IV}}(t) &= y^{\text{IV}}(t) \otimes \mathcal{L}^{-1}\{F_{\alpha_j}^{(\text{iter})}(s)\}, \end{aligned} \right. \quad (24)$$

with $y^{\text{IV}} = \frac{\hat{B}(p)}{A(p)}u(t)$. The regression vectors $\varphi_f^*(t)$ and $\varphi_f^{\text{IV}}(t)$ are formed as in (12) and (19). Finally, the parameter vector $\hat{\theta}_{\text{srivcf}}^{(\text{iter})}$ is computed at each iteration as:

$$\hat{\theta}_{\text{srivcf}}^{(\text{iter})} = (\mathbf{\Phi}_f^{\text{IV}^T} \mathbf{\Phi}_f^*)^{-1} \mathbf{\Phi}_f^{\text{IV}^T} \mathbf{Y}_f^*, \quad (25)$$

with

$$\begin{aligned} \mathbf{\Phi}_f^{\text{IV}^T} &= [\varphi_f^{\text{IV}}(T_s) \quad \varphi_f^{\text{IV}}(2T_s) \quad \dots \quad \varphi_f^{\text{IV}}(KT_s)]^T \\ \varphi_f^{\text{IV}}(t)^T &= \begin{bmatrix} \mathbf{D}^{\beta_0} u_f(t) \cdots \mathbf{D}^{\beta_{m_B}} u_f(t) \\ -\mathbf{D}^{\alpha_1} y_f^{\text{IV}}(t) \cdots -\mathbf{D}^{\alpha_{m_A}} y_f^{\text{IV}}(t) \end{bmatrix}. \end{aligned} \quad (26)$$

As in the rational case, the srivcf estimator is asymptotically unbiased whatever additive zero-mean noise is present. When the additive noise is white, the srivcf parametric error covariance matrix:

$$P_{\hat{\theta}} = \hat{\sigma}^2 (\mathbf{\Phi}_f^{\text{IV}^T} \mathbf{\Phi}_f^{\text{IV}})^{-1}, \quad (27)$$

is minimal, with $\hat{\sigma}^2$ the estimation of the noise variance.

### 2.2.2. Optimality of SRIVCF estimation
In any parameter estimation problem, the estimator performance can be assessed by comparing the estimator covariance matrix to the best possible accuracy given by the Cramér-Rao lower-bound (CRB) (see Kay (1993); Stoica and Moses (2005); Larsson et al. (2008)). Under a Gaussian noise assumption, which is the case here, the CRB is computable and often tight, meaning that there exists an IV estimator (Wellstead (1978)) that asymptotically reaches the CRB. In this regard, it has been shown in Söderström and Stoica (1989) that, when the noise model $H_o(p)$ is independent

of $\theta$, the minimum value of the covariance matrix $P_{\hat{\theta}}$, as a function of the design variables $\Phi_f^{\text{IV}}$ and $F(p)$ exists and is given by

$$P_{\hat{\theta}} \geq P_{\theta}^{\text{opt}}. \quad (28)$$

$P_{\theta}^{\text{opt}} = \sigma^2 \text{E} \left\{ \Phi_f^{\text{IV}} \Phi_f^{\text{IV}^T} \right\}$ is obtained with the true design variables. $\Phi_f^{\text{IV}}$ is the optimally noise-free pre-filtered IV vector, with the associated design variables: $F(p) = \frac{1}{H_o(p)A(p)}$, $H_o(p) = 1$ for a white Gaussian noise and $\Phi_f^{\text{IV}}(t) = F(p) [\mathbf{D}^{\beta_0}u(t) \cdots \mathbf{D}^{\beta_{m_B}}u(t) - \mathbf{D}^{\alpha_1}y^{\text{IV}}(t) \cdots - \mathbf{D}^{\alpha_{m_A}}y^{\text{IV}}(t)]^T$. The reader can refer to Söderström and Stoica (1989) for the proof of this bound.

### 2.2.3. Summary of the srivcf algorithm
The overall srivcf algorithm is summarized below:

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

**Step 2** iterative instrumental variable estimation with prefilters.

for iter = 2 to convergence

(i) Generate the instrumental variables $y^{\text{IV}}$ from the auxiliary model with the estimated polynomials based on the estimated parameter vector $\theta^{(\text{iter}-1)}$.

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

(iii) Based on these prefiltered data, compute the new estimates $\theta^{(\text{iter})}$ as in (25).

end for

**Step 3** compute the estimated parametric error covariance matrix associated with the parameter estimates with (27).

For a given $\gamma$, the parametric vector $\theta$ converges to an optimal value with an error criterion

$$J(\gamma) = \varepsilon^T \varepsilon, \quad (29)$$

with

$$\begin{aligned} \varepsilon &= y^*(t) - \hat{y}, \\ \hat{y} &= \frac{\hat{B}(p)}{\hat{A}(p)}u(t). \end{aligned} \quad (30)$$

With a fixed parameter vector, the next step consists in optimizing the commensurate order $\gamma$.

## 2.3 Commensurate order optimization

Several gradient-based methods exist in the literature (see Ljung (1999)), here the Gauss-Newton (GN) method is used. The commensurate order $\gamma_k$ at iteration $k$ depends on its previous value $\gamma_{k-1}$ and on a correction term:

$$\gamma_k = \gamma_{k-1} - \lambda H^{-1} \frac{\partial J}{\partial \gamma} \Big|_{\gamma_{k-1}}, \quad (31)$$

with $\lambda$ an evolutive factor that will accelerate the convergence to a minimum and
$$ \left. \frac{\partial J}{\partial \gamma} \right|_{\gamma_{k-1}} = 2 \left. \frac{\partial \varepsilon}{\partial \gamma} \right|^{T} \varepsilon, \quad H = 2 \left. \frac{\partial \varepsilon}{\partial \gamma} \frac{\partial \varepsilon}{\partial \gamma} \right|_{\gamma_{k-1}}, \quad \left. \frac{\partial \varepsilon}{\partial \gamma} \right| = - \frac{\partial \hat{y}}{\partial \gamma}, $$

$$ \frac{\partial \hat{y}}{\partial \gamma} = u(t) \otimes $$

$$ \mathcal{L}^{-1} \left\{ \ln(s) \frac{\sum_{i=0}^{n} i \hat{b}_i s^{i\gamma} + \sum_{i=0}^{n} \sum_{j=1}^{m} (i-j) \hat{b}_i \hat{a}_j s^{(i+j)\gamma}}{\left(1 + \sum_{j=1}^{n} \hat{a}_j s^{j\gamma}\right)^2} \right\}. $$

However $\frac{\partial \hat{y}}{\partial \gamma}$ is function of $\ln(s)$ which is quite difficult to simulate. Thus, the numerical derivative of $\frac{\partial \varepsilon}{\partial \gamma}$ is computed. Computing the approximated Hessian $H$, instead of the Hessian, is quicker as the sensitivity $\frac{\partial \hat{y}}{\partial \gamma}$ function is already computed for the gradient.

## 2.4 Summary of the Commensurate Order Optimization

The algorithm optimizes the commensurate order by using the GN method and the srivcf algorithm. The proposed algorithm enables to converge to a local minimum by minimizing the error criterion (29):

**Step 1** Initialize $\gamma_0$ and $\Lambda$

Compute $J(\gamma_0)$

$k = 1$

**Step 2** Gauss-Newton Optimization method
do

$$ \gamma_{temp} = \gamma_{k-1} - \lambda H^{-1} \left. \frac{\partial J}{\partial \gamma} \right|_{\gamma_{k-1}} $$

Compute the error criterion: $J(\gamma_{temp})$.

if $J(\gamma_{temp}) < J(\gamma_{k-1})$

$$ k = k + 1 $$

$$ \gamma_k = \gamma_{temp} $$

$$ \lambda = \Lambda $$

else

$$ \lambda = \lambda/2 $$

end if

$$ \text{while } |\gamma_k - \gamma_{k-1}| > \epsilon_1 \text{ and } |J(\gamma_k) - J(\gamma_{k-1})| > \epsilon_2 $$

Let us consider that the working area of the criterion $J$ is convex. At iteration $k$, it may happen that $J(\gamma_k) > J(\gamma_{k-1})$, but there exist a local minimum that dwells in $[\gamma_{k-1}, \gamma_k]$ (if $\gamma_{k-1} < \gamma_k$). Because of a $\lambda$ too high, the criterion $J(\gamma_k)$ is higher than $J(\gamma_{k-1})$. The criterion $J(\gamma_k)$ is computed again with a reduced convergence parameter $\lambda/2$ so that the algorithm converges to that local minimum when $J(\gamma_k) < J(\gamma_{k-1})$ (or until a precision on $J(\gamma_k)$ is obtained).

# 3. NUMERICAL EXAMPLE

## 3.1 System specifications

To illustrate the efficiency of this algorithm in system identification using fractional models, the following transfer function is used to generate data:

$$ G_0(s) = \frac{K(-Ts^{\gamma} + 1)}{\left(\left(\frac{s}{\omega_1}\right)^{2\gamma} + 2\zeta_1\left(\frac{s}{\omega_1}\right)^{\gamma} + 1\right)\left(\left(\frac{s}{\omega_2}\right)^{2\gamma} + 2\zeta_2\left(\frac{s}{\omega_2}\right)^{\gamma} + 1\right)}, \quad (32) $$

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

<figure><img src="image_4.png" /><figcaption>Fig. 3. Input and noisy output signals of the system model.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Error criteria computed with srivcf.</figcaption></figure>

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

The input signal $u(t)$, a pseudo random binary sequence of amplitude $[-5; 5]$, is applied. The output $y(t)$ is corrupted by an additive Gaussian white noise $p(t)$ with zero mean and a signal to noise ratio of 20dB as shown in Fig. 3. The sampling period is fixed to $T_s = 5.10^{-2}s$.

The instrumental variables are computed according to the iterative algorithm described in §2.2.3 in order to ripen the parametric estimation. The model, which parameters are to be estimated, is set to the same class as $G_0(s)$:

$$ G(s) = \frac{b_1 s^{\gamma} + b_0}{a_4 s^{4\gamma} + a_3 s^{3\gamma} + a_1 s^{2\gamma} + a_5 s^{\gamma} + 1}. \quad (33) $$

An optimization of the commensurate order is done in order to minimize the following normalized error criterion:

$$ J = \frac{\varepsilon^T \varepsilon}{y^T y}, \quad (34) $$
<figure><img src="image_6.png" /><figcaption>Fig. 5. (a) Histogram of the estimated commensurate order for 200 MCS runs with *coosrivcf* (b) Poles estimated for the 200 MCS where + denote the true poles</figcaption></figure>

with $\epsilon = y^* - \hat{y}$, and $\hat{y}$ is the estimated output. The Gauss-Newton method converges to the closest minimum. The *coosrivcf* algorithm, resulting from the approximated Hessian, optimizes the commensurate order and will reach the global minimum if the criterion *J* is convex. Fig. 4 illustrates the global error criterion: at each $\gamma$ varying in [0, 2] with a step of 0.05, the *srivcf* is applied and the error criterion is computed. The *coosrivcf* enables to obtain the global minimum unless working in a non convex area.

A Monte Carlo simulation (MCS) of 200 runs has been undertaken with $\gamma_0 = 0.3$. As expected, the *coosrivcf* algorithm converges to the closest minimum (here, true commensurate order 0.5). Fig. 5-a and Fig. 5-b present the accuracy of this algorithm as well for the estimation of the commensurate order and the consistency of the parameters by plotting the estimated poles. With a mean value of $\bar{\gamma} = 0.5004$ and a standard deviation of $\sigma(\gamma) = 0.013$, the MCS shows that the estimation is very accurate. Knowing the performance and efficiency of the *srivcf* (see Malti et al. (2008b) for comparisons with other methods) and being in the same model class as $\gamma = 0.5$, *srivcf* gives consistent estimates.

With an initial value $\gamma_0 = 0.3$ and a given noise realization, the *coosrivcf* algorithm converged to the following model:

$$ \hat{G}(s) = \frac{0.502s^{0.499} - 1.000}{24.999s^{1.99} + 36.521s^{1.498} + 31.178s^{0.999} + 5.311s^{0.499} + 1} \quad (35) $$

The estimated transfer function $\hat{G}(s)$ converges to the true system $G_0(s)$. Thus, the model is validated.

## 4. CONCLUSION

In this paper, the simplified refined instrumental variable for continuous-time systems for fractional models has been recalled (*srivcf*). This algorithm computes consistent estimates with minimum variance when being in the same class as the true system which requires the commensurate order to be known. If unknown, the commensurate order

has to be optimized, which is implemented in the *coosrivcf* algorithm. As shown in the example, the cost function is tightly linked to the commensurate order. The parametric estimation has been extended to the commensurate order.

## REFERENCES

M. Aoun, R. Malti, F. Levron, 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.

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

S.M. Kay. *Fundamentals of statistical signal processing: estimation theory*. Prentice-Hall, Upper Saddle River, NJ, 1993.

E.K. Larsson, M. Mossberg, and T. Södreström. *Identification of continuous-time models from sampled data*, chapter Estimation of continuous-time Stochastic System Parameters. Springer-Verlag, H. Garnier and L. Wang (Eds.) edition, 2008.

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

R. Malti, X. Moreau, and F. Khemane. Resonance of fractional transfer functions of the second kind. In *3rd IFAC workshop on fractional differentiation and its applications FDA'08*, Ankara, Turkey, 2008a.

R. Malti, S. Victor, A. Oustaloup, and H. Garnier. An optimal instrumental variable method for continuous-time fractional model identification. In *17th World Congress IFAC'08*, Seoul, Korea, July 2008b.

K.B. Oldham and J. Spanier. *The fractional calculus. Academic Press, New-York and London, 1974.*

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

T. Poinot and J. C. Trigeassou. A method for modelling and simulation of fractional systems. *Signal processing*, 83:2319-2333, 2003.

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*. Prentice Hall, 1989.

P. Stoica and R. Moses. *Spectral Analysis of Aignals*. Pearson Prentice-Hall, Upper Saddle River, NJ, 2005.

P. E. Wellstead. An instrumental product moment test for model order estimation. *Automatica*, 14:89-91, 1978.

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

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

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.