# From System Identification to Path Planning Using Fractional Approach: A Thermal Application Example

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

► To cite this version:

Stéphane Victor, Rachid Malti, Pierre Melchior, Alain Oustaloup. From System Identification to Path Planning Using Fractional Approach: A Thermal Application Example. ASME 2009 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Aug 2009, San Diego, United States. pp.235-242, 10.1115/DETC2009-87014. hal-05447134

HAL Id: hal-05447134

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

Submitted on 10 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
FROM THE IDENTIFICATION TO THE TRAJECTORY PLANNING OF A THERMAL
APPLICATION

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

IMS, UMR 5218 CNRS, Université Bordeaux, 351 cours de la Libération,
F 33405 Talence cedex, France –

Email: {stephane.victor, rachid.malti, alain.oustaloup}@laps.ims-bordeaux.fr

ABSTRACT

Given a specific trajectory, the objective of this paper is to
follow a desired path using fractional linear flatness. However,
models of true system are not always known. Through a black
box fractional identification, the heat flux of an aluminium rod
can be modeled. The coosrivcf enables to optimize the commen-
surate order and the fractional flat output generates the input
command so that the system follows a desired path.

# 1 Introduction

Theoretical and experimental studies point out that thermal
[1], electro-chemical [2, 3] and viscoelastic [4] systems are de-
fined by non integer differential equations. Integer models are
not appropriated to identify that type of systems, unless apply-
ing very complex models. As detailed in the survey paper by [5]
and all references therein, **Instrumental Variable (IV)** approaches
to system identification and parameter estimation have been de-
veloped since 1960s and 1970s. Recently, **Refined Instrumental
Variables for Continuous-system** [6, 7] has been extended for
identifying fractional continuous-time systems [8]. [9] are the
first to extend the use of state variable filters to fractional sys-
tems, with an IV procedure to eliminate bias. However, their
instruments are not selected in an optimal way. The optimal ap-
proach, *srivcf*, involves a method of iterative-adaptive prefilter-
ing based on a quasi-optimal statistical solution obtained from a
given commensurable<sup>1</sup> order. This optimal IV method presents

the advantage of yielding asymptotically consistent estimates in
noisy context, 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 commensurable differentiation
order. These are however usually unknown in practice which
might induce unmodeled dynamics. Consequently, after fixing
the dimension of the system, it is proposed to combine the *srivcf*
algorithm to the optimization of the commensurate order, which
is a nonlinear gradient-based method: the Levenberg-Marquardt
approach.

A new property called flatness characterizing a certain class
of nonlinear systems, allowing the achievement of a simple and
robust control, has been introduced by M. Fliess, J. Levine, Ph.
Martin and P. Rouchon (see [10, 11]). If we restrict to linear
time-invariant systems (LTI), it turns out that such systems are
flat if and only if they are controllable and that a particular flat
output may be obtained as a by-product of the Brunowsky con-
trollability canonical form. Flat outputs summarize, on its own,
all the system dynamics. All system variables will be deduced
from them without integrating differential equations. In previous
works on linear fractional flatness, the extension of the flatness
principle to fractional linear systems was established in [12] with
a transfer function approach, [13, 14] with a polynomial matrix
approach and state-space-like representation where the system
was supposed commensurate with a known order. In this paper,
the study involves trajectory planning of a system whose com-
mensurate order is unknown. Flatness is well-suited to trajectory
planning and offers the possibility to solve straight-forwardly
path tracking. One of the interest of flatness relies on computing

<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>
directly the command trajectory enabling to generate the desired trajectory without disturbance.

The objective here is to investigate thermal diffusion in an aluminium rod. From experimental data, a performant fractional model is determined using black box identification which estimates as well the coefficient and the commensurate order. Then, through fractional flatness principles, a flat output will enable to follow a desired path.

## 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$ [15]:

$$\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 negatives 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, see [16]:

$$\mathcal{L}\{\mathbf{D}^{\gamma} x(t)\} = s^{\gamma} X(s) \quad \text{if } x(t) = 0 \quad \forall t \le 0. \quad (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}}. \quad (5)$$

Moreover if $F(s)$ is commensurable 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}}, \quad (6)$$

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

$$\begin{cases} \tilde{b}_i = b_i & \text{if } i\gamma = \beta_i \quad \text{and} \quad \tilde{b}_i = 0 & \text{if } i\gamma \neq \beta_i \\ \tilde{a}_j = a_j & \text{if } j\gamma = \alpha_j \quad \text{and} \quad \tilde{a}_j = 0 & \text{if } j\gamma \neq \alpha_j. \end{cases} \quad (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 commensurable order $\gamma_k$ which is estimated iteratively.

By considering a first level of generalization, the derivation orders of each elementary differential equation are the same and given by the commensurate order, which is a scalar, such as:

$$\left\{ \begin{aligned} x^{(\nu)}(t) &= Ax(t) + Bu(t) \\ y(t) &= Cx(t) + Du(t) \end{aligned} \right. , \quad (8)$$

where $A, B, C$ and $D$ are respectively a $n \times n$, $n \times m$, $r \times n$ and $r \times m$ constant matrices.

## 2 Identification by fractional model

The optimization of the commensurable order is combined with the srivcf algorithm in a two stage procedure. First a commensurable order $\gamma_0$ is fixed. Then, for a fixed commensurable order, optimal coefficients are computed iteratively using the
<figure><img src="image_2.png" /><figcaption>Figure 1. Iterative order optimization algorithm: coosrivcf</figcaption></figure>

srivcf algorithm. Next, for the computed coefficients, an optimal commensurable order is calculated using nonlinear programming technics. This procedure is repeated until convergence as sketched in Fig. 1. Although such algorithms have no proof of convergence, it is noted that they do often converge.

The Refined Instrumental Variables for Continuous system [17] (rivc) method and its simplified version srivc (when the additive measurement noise is assumed to be white) can be interpreted in optimal statistical terms 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.

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

This stochastic approach proposed by [5, 6], srivc, includes an adaptive pre-filtering based on the optimal solution when the noise $\xi$ is a white gaussian one. This method is a logical extension of the instrumental variable estimator with SVF. This estimator is extended to fractional models and is named srivcf (srivc for fractional models) and where the commensurate order is known. Following the Prediction Error Minimization (PEM) approach, a suitable error function $\epsilon(t)$ in CT is given by the output error (OE) (by introducing the differential operator $p = \frac{d}{dt}$):

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

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

$$ \begin{aligned} \epsilon(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 (10) $$

As shown by [5], the optimal estimator (asymptotically without bias and of minimal variance) is obtained when the SVF filters $F_{\gamma}(s)$ are replaced by:

$$ F_{\gamma}^{\text{opt}}(p) = \frac{p^{\gamma}}{H(p)A(p)}, \quad (11) $$

where $H(p)$ is the real noise model, $A(p)$ 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.

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:

$$ \epsilon(t) = A(p)y_f(t) - B(p)u_f(t) \quad (12) $$

or,

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

where

$$ \left\{ \begin{array}{l} \mathbf{D}^{\beta_i} u_f(t) = u(t) * \mathcal{L}^{-1}\{F_{\beta_i}(s)\}, \quad i = 0, \dots, m_B \\ \mathbf{D}^{\alpha_j} y_f(t) = y^*(t) * \mathcal{L}^{-1}\{F_{\alpha_j}(s)\}, \quad j = 1, \dots, m_A \end{array} \right. \quad (14) $$

In practical cases, $A(s)$ is not known. The filter $F(s)$ will then be updated by an iterative procedure to estimate the denominator $\hat{A}(s)$. This algorithm will enable to optimize the instruments. Initialized by the $\hat{\theta}_{IVSVF}$ estimates (see [8] 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_{\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 (15) $$
where (iter) = 1, 2, ... stands for the iteration number and $\hat{a}_j$ the estimated $a_j$ coefficient at iteration (iter).

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) * \mathcal{L}^{-1}\{F_{\beta_i}^{(\text{iter})}(s)\} \\
\mathbf{D}^{\alpha_j} y_f^*(t) &= y^*(t) * \mathcal{L}^{-1}\{F_{\alpha_j}^{(\text{iter})}(s)\} \\
\mathbf{D}^{\alpha_j} x_f^{IV}(t) &= x^{IV}(t) * \mathcal{L}^{-1}\{F_{\alpha_j}^{(\text{iter})}(s)\},
\end{aligned}
\right. 
\quad (16) $$

where $x^{IV}(t) = \frac{B(s, \theta_{IV}^{(iter)})}{A(s, \theta_{IV}^{(iter)})}u(t)$. The regression vectors $\varphi_f^*(t)$ and $\varphi_f^{IV}(t)$ are formed as:

$$ 
\begin{aligned}
\varphi_f^*(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} \\
\varphi_f^{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} x_f^{IV}(t) \cdots -\mathbf{D}^{\alpha_{m_A}} x_f^{IV}(t) \end{bmatrix}.
\end{aligned} 
\quad (17) $$

Finally, the parameter vector $\hat{\theta}_{srivcf}^{(iter)}$ is computed at each iteration as:

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

with

$$ 
\begin{aligned}
Y_f^* &= [y_f^*(T_s) \quad y_f^*(2T_s) \quad \dots y_f^*(KT_s)]^T, \\
\Phi_f^* &= [\varphi_f^*(T_s) \quad \varphi_f^*(2T_s) \quad \dots \varphi_f^*(KT_s)]^T, \\
\Phi_f^{IV} &= [\varphi_f^{IV}(T_s) \quad \varphi_f^{IV}(2T_s) \quad \dots \varphi_f^{IV}(KT_s)]^T.
\end{aligned} 
\quad (19) $$

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 with the estimate of $\hat{\theta}_{srivcf}$ is given by:

$$ \mathbf{P}_{\hat{\theta}} = \hat{\sigma}^2 (\Phi_f^{IV T} \Phi_f^*)^{-1}, \quad (20) $$

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

At the $k^{th}$ sampling instant $t = t_k$, the associated estimation model can be written as

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

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

$$ \varphi_f(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 (22) $$

The estimation model (21) forms a basis for the definition of a likelihood function and maximum likelihood estimation.

## 2.2 Commensurable order optimization

Several gradient-based methods exist in the literature (see [18]), here the Levenberg-Marquardt (LM) method is used. The commensurable 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 (23) $$

$$ 
\begin{aligned}
J &= \varepsilon^T \varepsilon, \\
\left. \frac{\partial J}{\partial \gamma} \right|_{\gamma_{k-1}} &= 2 \left. \frac{\partial \varepsilon}{\partial \gamma} \right)^T \varepsilon, \\
H &= 2 \left. \frac{\partial \varepsilon}{\partial \gamma} \frac{\partial \varepsilon}{\partial \gamma} \right|_{\gamma_{k-1}}, \\
\frac{\partial \varepsilon}{\partial \gamma} &= - \frac{\partial \hat{\gamma}}{\partial \gamma}, \\
\frac{\partial \hat{\gamma}}{\partial \gamma} &= u(t) \otimes \mathcal{L}^{-1} \left\{ \ln(s) \frac{\sum_{i=0}^{n} \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\}.
\end{aligned}
 $$

However $\frac{\partial \hat{\gamma}}{\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. The advantage is on the computation time of the approximated Hessian $H$, as the second derivative of $\varepsilon$ is not necessary and avoids numerical problems.

## 2.3 Summary of the Commensurable Order Optimization

The proposed algorithm enables to converge to a local minimum by minimizing the error criterion $J$. Give an initial $\Lambda$ and a commensurable order $\gamma_0$ to initialize the algorithm:

**Step 1** Initialize $\gamma_0$
Compute $J(\gamma_0)$
Set $k = 1$

**Step 2** Levenberg-Marquardt Optimization method

**do**

Initialize $\lambda = \Lambda$

**do**

(i) Affine the order to converge to a minimum:

$$ \gamma_k = \gamma_{k-1} - \lambda H^{-1} \frac{\partial J}{\partial \gamma} \Big|_{\gamma_{k-1}} $$
(ii) Compute the error criterion: $J(\gamma_k)$.

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

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

$k = k + 1$

while $|\gamma_k - \gamma_{k-1}| > \epsilon_2$

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 Fractional flat output from polynomial forms

All along this paper, $u^{(k)}$ stands for the $k$-th order derivative of $u$ with respect to time, $k$ being a non integer scalar and each component of $u$ is derived at order $k$. Let the system be given by the following differential equation:

$$x^{(v)} = f(x, u), \quad (24)$$

with $x \in \mathbb{R}^n$, $u \in \mathbb{R}^m$ being the input command, $v \in [0, 2]^2$ and $f = (f_1, \dots, f_n)$ a regular function of $x$ and $u$ where the rank of $\frac{\partial^v f}{\partial u^v}$ is equal to $m$.

This system is said differentially flat (see [13,20]) if one can find a set of independent variables referred to as the flat output $z$ such that each system variable, including the inputs, is a function of $z$ and a finite number of its time derivatives. The fractional flat output $z$ is such that

$$z = h(x, u, u^{(v)}, u^{(2v)}, \dots, u^{(\beta v)}), \quad z \in \mathbb{R}^m, \quad (25)$$

with $\beta$ a finite $m$-tuple of integers $^3$ and $v$ stands for the commensurable order.

#### 3.1 Fractional flat output of a fractional linear system

In the following, the Laplace transformation $\mathcal{L}$, which is an algebraic tool generally used to represent fractional differentiation systems, is defined by (see [16]):

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

All variables are in the operational domain; in order to get back to the time domain, the Laplacian inverse operator $\mathcal{L}^{-1}$ should

be applied. Fractional linear systems such as (8) can be rewritten in the following polynomial matrix representation:

$$\left\{ \begin{array}{l} A(s^v)x = Bu \\ y = Cx + Du \end{array} \right. , \quad (27)$$

with $A(s^v) = s^v I - A$ an $n \times n$ matrix whose entries are polynomials of the formal variable $s$ ($I$ being the identity matrix), and $B, C$ and $D$ are constant matrices. System (27) is assumed controllable, i.e. $A$ and $B$ are left coprime (see [11]).

A fractional linear flat output defines an output defined by (25) with $h$ being linear; the state variables, the input command and the output can be redefined with $z$ such as:

$$\begin{align} x_i &= \sum_{j=1}^{m} \sum_{k=0}^{\alpha_j} a_{i,j,k} z_j^{(kv)}, & i &= 1, \dots, n, \notag \\ u_l &= \sum_{j=1}^{m} \sum_{k=0}^{\alpha_{j+1}} b_{l,j,k} z_j^{(kv)}, & l &= 1, \dots, m, \tag{28} \\ y_q &= \sum_{j=1}^{m} \sum_{k=0}^{\sigma_j} c_{q,j,k} z_j^{(kv)}, & q &= 1, \dots, r, \notag \end{align}$$

where $\alpha$ and $\sigma$ are $m$-tuple of integers. Notice that all differential orders are multiple of the commensurable order $v$.

In the polynomial matrix language, (28) reads

$$x = P(s^v) z, \quad u = Q(s^v) z, \quad (29)$$

where $P$ (resp. $Q$) is an $n \times m$ (resp. $m \times m$) polynomial matrix of $s^v$, with entries $P_{i,j}(s^v) = \sum_{k=0}^{q_j} a_{i,j,k} s^{v_{i,k}}$ (resp. $Q_{l,j}(s^v) = \sum_{k=0}^{q_j+1} b_{l,j,k} s^{v_{l,k}}$). For an overview of $s^v$-polynomials algebra, refer to [13].

Matrices $P$ and $Q$ satisfying (29) are called defining matrices of the fractional flat output $z$.

**Theorem 3.1.** The variable $z = (z_1, \dots, z_m)$ is a linear flat output of system (27) if and only if its defining matrices $P$ and $Q$ are given by

$$\begin{align} R^T A(s^v) P(s^v) &= 0, \tag{30} \\ A(s^v) P(s^v) &= BQ(s^v), \notag \end{align}$$

with $R$ an arbitrary matrix of rank $n-m$ orthogonal to $B$ (it means that $R^T B = 0$), and with $P(s^v)$ and $Q(s^v)$ of rank $m$ for every $s^v$ and right coprime $^4$.

In addition, linear flat output $y$ of the controllable system (27) always exists (and therefore its defining matrices $P$ and $Q$ always exist too).

<footnote>$^2$It is proved in [19] that fractional differentiation systems are unstable for a commensurate order $v \notin [0, 2]$.</footnote>

<footnote>$^3$For a $m$-tuple of integers $\beta = (\beta_1, \dots, \beta_m)$ of dimension $m$ and a vector $u$ of dimension $m$, the notation $u^{(\beta)}$ represents $(u^{(\beta_1)}, \dots, u^{(\beta_m)})$.</footnote>

<footnote>$^4$For the proof for fractional systems, refer to [13].</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2. Insulated long aluminium rod heated by a resistor</figcaption></figure>

## 3.2 Dynamic inversion of a fractional system

Once the state variable $x$ (resp. the input command $u$) being expressed by the flat output $z$ with its defining matrix $P$ (resp. $Q$), the desired output $y_{des}$ can be expressed with $z$. The fractional differentiation system is inverted in order to get the fractional flat output. From (8) and (29), the desired output yields:

$$y_{des} = (CP(s^v) + DQ(s^v))z. \quad (31)$$

Thus, it is necessary to invert (31) in order to determine the flat output entirely. We note $W(s^v) = CP(s^v) + DQ(s^v)$ which is of dimension $r \times m$.

1) if $r = m$ and rank($W$) = $m$, $W$ is reversible, and the flat output is easily determined as: $z = W^{-1}y_{des}$.

2) if $r < m$, $W$ is left reversible, and the flat output can be expressed such as: $z = (W^T W)^{-1} W^T y_{des}$.

3) if $r > m$, in this case, the system has more outputs than input commands. Some outputs are linear combination of independent outputs or linear combination of the command inputs. After isolating those outputs, the system can be brought to one of the two previous cases.

# 4 Physical modeling of the aluminium rod

## 4.1 Plant description

A long aluminium rod heated by a resistor is considered in this experiment. To ensure unidirectional heat transfer, the entire surface of the rod is insulated by a foam as shown in Fig. 2. The input signal is a thermal flux generated by a resistor glued at one end and the output signal is the temperature of the rod measured at a distance $x^* = 0.5cm$ from the heated end. The injected heat flux is controlled by a computer through an on-off transistor with a controlled amplitude of the input voltage.

## 4.2 Physical modeling

To demonstrate the fractional behavior of this thermal system, the rod is modeled under the following assumptions:

(i) the rod is perfectly isolated,

(ii) the rod is considered as a semi-infinite homogeneous plane medium with conductivity $\lambda$ and diffusivity $\alpha$,

(iii) at rest, the rod is at ambient temperature, so that there is no thermal exchange with the environment,

(iv) losses on the surface where the thermal flux is applied are neglected. The electrical energy consumed in the resistor is assumed to be totally transformed into thermal energy and diffused by conduction in unidirectionally.

Care is taken so that assumption (iii) is fulfilled. A time period of 24 hours is observed between two consecutive experiments so that the temperature of the aluminium rod cools down to the ambient temperature and hence initial conditions are zero at the beginning of each experiment, since there is initially no thermal exchange with the environment. Fractional systems with non zero initial conditions are much harder to handle (see [21, 22]).

The injected thermal flux is a pseudo random binary sequence (PRBS) varying from 0 to $\phi \approx 41kWm^{-2}$. The interest of using fractional models is to model systems that have fractional behaviors by using few parameters. It has been shown in [8, 23] that high order rational models are necessary to get comparable results to fractional models. Let us assume that the commensurate order is not known.

## 4.3 Black box identification

The equation error model is now used for the identification of the aluminium rod based on measurement data. The dimension of the system is fixed, as is usually done in the rational case. The procedure as sketched in Fig. 1 is applied. A time lag is noticed between the input and the output signals, probably generated by the flux diffusion in the medium. Consequently, Based on the study of the physical model, the following model structure is chosen:

$$F_1(s) = \frac{b_0}{a_3 s^{3\gamma} + a_2 s^{2\gamma} + s^{\gamma}} \times e^{-3s}. \quad (32)$$

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

$$J_{dB} = 10 \log_{10} \frac{\varepsilon^T \varepsilon}{y^T y}, \quad (33)$$

with $\varepsilon = y^* - \hat{y}$, and $\hat{y}$ is the estimated output. The coosrivcf method, resulting from the approximated Hessian, optimizes the commensurate order and will reach the global minimum if the criterion $J$ is convex. The coosrivcf enables to obtain the global minimum unless working in a non convex area.
<figure><img src="image_4.png" /><figcaption>Figure 3. Measured output compared to the fractional model output on identification data.</figcaption></figure>

With an initial $γ_0 = 0.1$, a $Λ = 10$, coosrivcf converged to the following model:

$$G(s) = \frac{1.394 \cdot 10^{-5}}{s^{3\gamma} - 0.069s^{2\gamma} + 0.058s^{\gamma}} \times e^{-3s}, \quad (34)$$

where $γ = 0.26$. To show the efficiency of the coosrivcf algorithm, the identified output and the non-noised measured output are plotted in Fig. 3. The identified model fit exactly to the experimental data.

## 4.4 Linear fractional flatness

The aim here is to control the temperature at a specific point of the rod $x^* = 0.05m$. With (34), the pseudo-state-space representation is as follows:

$$A = \begin{bmatrix} 0.069 & -0.058 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}, C = \begin{bmatrix} 0 & 0 & 1.394 \cdot 10^{-5} \end{bmatrix} \quad (35)$$

$$B = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}^T, D = 0$$

Let us apply the fractional flatness principles to trajectory planning. It is assumed that the first and second derivatives of the temperature are equal to zero at initial and final temperatures. The metallic rod will have to rise 30°C above the initial temperature in 1250 seconds and will stabilize for the same period of time, according to the following chosen path, a Polynomial Interpolation of degree (PI5):

$$y_n(t) = q_i + 80(q_{f,n} - q_i) \frac{t^3}{t_f^3} - 240(q_{f,n} - q_i) \frac{t^4}{t_f^4} + 192(q_f - q_i) \frac{t^5}{t_f^5}, \quad (36)$$

<figure><img src="image_5.png" /><figcaption>Figure 4. Simulation of the metallic rod in open-loop: input commands ($u = \Phi$), flat output, desired and effective (dashed) outputs</figcaption></figure>

with $q_i = 0°C$, $q_{f,n=1} = 30°C$, and $t_f = 2500s$. The desired output is illustrated in full line in the third plot of Fig. 4. (35) has the form (27) with

$$A(s) = \begin{bmatrix} s^\nu - 0.069 & 0.058 & 0 \\ -1 & s^\nu & 0 \\ 0 & -1 & s^\nu \end{bmatrix}. \quad (37)$$

Let us compute a matrix $R$ orthogonal to $B$:

$$R^T = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad (38)$$

$$R^T A(s) = \begin{bmatrix} -1 & s^\nu & 0 \\ 0 & -1 & s^\nu \end{bmatrix}. \quad (39)$$

Thus, according to Theorem 3.1, a flat output is obtained by computing $P(s)$ solution to $R^T A(s) P(s) = \begin{bmatrix} -P_1 + s^\nu P_2 \\ -P_2 + s^\nu P_3 \end{bmatrix} = 0$, which yields $P(s) = [s^{2\nu} s^\nu 1]^T$. $P(s)$ is clearly of rank 1 for all $s^\nu$. Finally, the input command $u$ (to be applied to the system in order to get the desired output), and the state $x$ are computed according to (29). Thus, the effective output is also generated with (31). Fig. 4 illustrates the input command $u$, the fractional flat output (which has no physical meaning) and the effective output which is compared to the desired path.

## 5 Conclusion

The optimal instrumental variable estimator, srivcf, is recalled as it has the advantage to yield asymptotically consistent
estimates and is of minimum variance in presence of white noise. Moreover, in a black box identification, the coosrivcf, algorithm has enabled to estimate as well the model parameters and the commensurate order. These were applied to real experimental data: the heat transfer in an aluminium rod. The fractional identification is well suited here as it requires less parameters to achieve a good quality of identification. Moreover, through fractional flatness principles, a desired trajectory could have been planned to the identified fractional model. The fractional flat output has enabled to generate the input command that yield the desired temperature path. Further studies could be done in a closed-loop approach with robust trajectory planning in order to reject perturbations or with MIMO applications.

## REFERENCES

[1] Battaglia, J.-L., Le Lay, L., Batsale, J.-C., Oustaloup, A., and Cois, O., 2000. "Heat flux estimation through inverted non integer identification models". *International Journal of Thermal Science*, **39**(3), pp. 374-389.

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

[3] Sabatier, J., Aoun, M., Oustaloup, A., Grégoire, G., Ragot, F., and Roy, P., 2006. "Fractional system identification for lead acid battery state charge estimation". pp. 2645-2657.

[4] Moreau, X., Ramus-Serment, C., and Oustaloup, A., 2002. "Fractional differentiation in passive vibration control". *International Journal of Nonlinear Dynamics and Chaos in Engineering Systems*, **Vol. 29**(1-4), July, pp. 343-362. Special Issue on "Fractional Order Calculus and its Applications".

[5] Young, P., 1981. "Parameter estimation for continuous-time models - a survey". *Automatica*, **17**(1), pp. 23-29.

[6] Young, P., 2002. "Optimal IV identification and estimation of continuous-time TF models". In 15th World IFAC Congress, IFAC, Elsevier.

[7] Garnier, H., 2006. "Identification de modèles paramétriques à temps continu. méthodes, outli logiciel et avantages". *Journées Identification et Modélisation Expérimentale (JIME'2006)*, november. Poitiers (France).

[8] Malti, R., Victor, S., and Oustaloup, A., 2008. "An optimal instrumental variable method for continuous-time fractional model identification". In 17th World Congress IFAC'08.

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

[10] Fliess, M., Lévine, J., Martin, P., and Rouchon, P., 1999. "A Lie-Bäcklund approach to equivalence and flatness of

nonlinear systems". *IEEE Trans. of Autom. Control*, **44**, pp. 922-937.

[11] Levine, J., and Nguyen, D. V., 2003. "Flat output characterization for linear systems using polynomial matrices". *System and Control Letters*, **48**, pp. 69-75.

[12] Melchior, P., Cugnet, M., Sabatier, J., Poty, A., and Oustaloup, A., Advances in Fractional Calculus Theoretical Developments and Applications in Physics and Engineering, Springer 2007. *Flatness control of a fractional thermal system*. In Sabatier, J. and Agrawal, O. P. and Tenreiro Machado, J. A. (Eds), ch. VII, pp. 493-509.

[13] Victor, S., Melchior, P., Sabatier, J., and Oustaloup, A., 2008. "Extension de la platitude aux systèmes fractionnaires MIMO : application un système thermique". In IEEE CIFA'08.

[14] Victor, S., Melchior, P., Nelason-Gruel, D., and Oustaloup, A., 2008. "Flatness control for linear fractional MIMO systems: thermal application". In 3rd IFAC Workshop on FDA'08.

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

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

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

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

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

[20] Fliess, M., Lévine, J., Martin, P., and Rouchon, P., 1995. "Flatness and defect of nonlinear systems: introductory theory and examples". *International Journal of Control*, **61**, pp. 1327-1361.

[21] Lorenzo, C., and Hartley, T., 2008. "Initialization of fractional-order operators and fractional differential equations". *Journal of computational and nonlinear dynamics. Special issue on discontinuous and fractional dynamical systems*, **3**, April, pp. 021101-1,021101-9.

[22] Sabatier, J., Merveillaut, M., Malti, R., and Oustaloup, A., 2008. "On a representation of fractional order systems: Interests for the initial condition problem". In 3rd IFAC workshop on FDA'08.

[23] Malti, R., Sabatier, J., and Hüseyin, A., 2009. "Thermal modeling and identification of an aluminium rod using fractional calculus". In submitted, 15th IFAC Symposium on System Identification (SYSID'2009).