Identification and uncertainty domain quantification: experimental validation on a "smart" wind turbine
Available from
Guillaume Mercère's profile on Mendeley.
Page 1
Identification and uncertainty domain quantification: experimental validation on a "smart" wind turbine
Identification and Uncertainty Domain
Quantification: Experimental Validation on
a ’Smart’ Wind Turbine
W. Farah ∗ G. Merce`re ∗ I. Houtzager ∗∗
J. W. Van Wingerden ∗∗ T. Poinot ∗
∗University of Poitiers, Laboratoire d’Automatique et Informatique
Industrielle, 2 rue Pierre Brousse, Baˆtiment B25, 86022 Poitiers,
France (e-mail: wafa.farah@etu.univ-poitiers.fr,
guillaume.mercere@univ-poitiers.fr, thierry.poinot@univ-poitiers.fr).
∗∗University of Delft, Delft Center of System and Control, Delft, 2628
CD The Netherlands, (e-mail: I.Houtzager@tudelft.nl,
J.W.vanWingerden@tudelft.nl)
Abstract: In this paper, a recent subspace-based identification approach is used to obtain an
experimental model of a smart rotor, which is a wind turbine rotor equipped with a number
of additional actuators and sensors to react to turbulence in a more detailed way. To make
sure that the eventual controller is robust, it is essential to quantify the model uncertainty. In
this paper a novel bounded error approach is tailored towards the characteristics of the smart
rotor. It is shown that the suggested approach can provide realistic uncertainty bounds on the
estimated model of a prototyped smart rotor.
Keywords: Parameter identification; uncertainty quantification; smart rotor.
1. INTRODUCTION
Identification for robust control must deliver not only a
nominal model, but also a reliable estimate of the un-
certainties associated with the model (Ljung, 1999a). To
reach this goal, two steps are generally necessary. The first
one deals with the identification of the model. This model
is required to understand, to control or to improve the
system’s functioning. The second step consists in design-
ing the control law from the estimated model. However,
because the model is only an approximation of the system’s
dynamics, it is paramount to fix some constraints so that
the controller designed from the identified model achieves
good performance when it is applied to the real system.
In other words, the controller must be robust with respect
to the uncertainties of the estimated parameters. Thus,
the uncertainties of the estimated model must be well-
described.
In system identification theory, the uncertainty domain
description is mainly based on prior assumptions about
noise and unmodeled dynamics (Goodwin et al., 1992;
Ljung, 1999a). In (Ljung, 1997, 1999b), it is assumed that
the disturbances acting on the system are realizations
of random variables. Because the information on the
measurement noise is seldom available and/or difficult to
be verified (de Vries, 1994), a different characterization
way can be considered (Ninness and Goodwin, 1995). This
approach is mainly based on deterministic hypotheses, i.e.,
on the assumption that the residuals are unknown-but-
bounded (Giarre´ and Milanese, 1997; Giarre´ et al., 1997;
Farah et al., 2010a,b). This basic idea has given rise to a
number of techniques usually addressed as bounded error
or set membership identification (Walter et al., 1996). The
main drawback of this approach is its dependence on the
way the bound is determined. Notice indeed that the error
comes from two different sources (the unmodeled dynamics
and the noise affecting the data) which makes the bound
determination quite difficult. To tackle this problem, some
efforts are made for estimating this modeling error (Ljung,
1999b; Reinelt et al., 2001).
In this paper, state-space model representations are used
in order to deal with SISO and MIMO systems in a similar
way. It is also well-known that such representations are
much more convenient for a future controller design. As
far as the identification problem is concerned, a modified
predictor-based subspace identification (PBSID) (Chiuso,
2007) adapted with the propagator method (Munier and
Delisle, 1991; Merce`re and Lovera, 2007; Merce`re et al.,
2008; Merce`re and Bako, 2011), is used in order to estimate
directly a state-space realization of the system from the
measured input-output (I/O) data. Contrary to the classic
subspace algorithms (Van Overschee and De Moor, 1996;
Katayama, 2005; Verhaegen and Verdult, 2007), this tech-
nique does not give access to a fully-parametrized form but
leads to a state-space representation with a minimal num-
ber of parameters, even for MIMO systems. In addition,
thanks to the propagator method, the state vector can be
estimated directly from the available I/O data (Houtzager
et al., 2009). By doing so, a model linear with respect to
the parameters (LP) is obtained which makes the descrip-
tion of uncertainty areas easier. These domains are derived
from the analysis of particular quadratic criteria.
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
Copyright by the
International Federation of Automatic Control (IFAC)
4404
Quantification: Experimental Validation on
a ’Smart’ Wind Turbine
W. Farah ∗ G. Merce`re ∗ I. Houtzager ∗∗
J. W. Van Wingerden ∗∗ T. Poinot ∗
∗University of Poitiers, Laboratoire d’Automatique et Informatique
Industrielle, 2 rue Pierre Brousse, Baˆtiment B25, 86022 Poitiers,
France (e-mail: wafa.farah@etu.univ-poitiers.fr,
guillaume.mercere@univ-poitiers.fr, thierry.poinot@univ-poitiers.fr).
∗∗University of Delft, Delft Center of System and Control, Delft, 2628
CD The Netherlands, (e-mail: I.Houtzager@tudelft.nl,
J.W.vanWingerden@tudelft.nl)
Abstract: In this paper, a recent subspace-based identification approach is used to obtain an
experimental model of a smart rotor, which is a wind turbine rotor equipped with a number
of additional actuators and sensors to react to turbulence in a more detailed way. To make
sure that the eventual controller is robust, it is essential to quantify the model uncertainty. In
this paper a novel bounded error approach is tailored towards the characteristics of the smart
rotor. It is shown that the suggested approach can provide realistic uncertainty bounds on the
estimated model of a prototyped smart rotor.
Keywords: Parameter identification; uncertainty quantification; smart rotor.
1. INTRODUCTION
Identification for robust control must deliver not only a
nominal model, but also a reliable estimate of the un-
certainties associated with the model (Ljung, 1999a). To
reach this goal, two steps are generally necessary. The first
one deals with the identification of the model. This model
is required to understand, to control or to improve the
system’s functioning. The second step consists in design-
ing the control law from the estimated model. However,
because the model is only an approximation of the system’s
dynamics, it is paramount to fix some constraints so that
the controller designed from the identified model achieves
good performance when it is applied to the real system.
In other words, the controller must be robust with respect
to the uncertainties of the estimated parameters. Thus,
the uncertainties of the estimated model must be well-
described.
In system identification theory, the uncertainty domain
description is mainly based on prior assumptions about
noise and unmodeled dynamics (Goodwin et al., 1992;
Ljung, 1999a). In (Ljung, 1997, 1999b), it is assumed that
the disturbances acting on the system are realizations
of random variables. Because the information on the
measurement noise is seldom available and/or difficult to
be verified (de Vries, 1994), a different characterization
way can be considered (Ninness and Goodwin, 1995). This
approach is mainly based on deterministic hypotheses, i.e.,
on the assumption that the residuals are unknown-but-
bounded (Giarre´ and Milanese, 1997; Giarre´ et al., 1997;
Farah et al., 2010a,b). This basic idea has given rise to a
number of techniques usually addressed as bounded error
or set membership identification (Walter et al., 1996). The
main drawback of this approach is its dependence on the
way the bound is determined. Notice indeed that the error
comes from two different sources (the unmodeled dynamics
and the noise affecting the data) which makes the bound
determination quite difficult. To tackle this problem, some
efforts are made for estimating this modeling error (Ljung,
1999b; Reinelt et al., 2001).
In this paper, state-space model representations are used
in order to deal with SISO and MIMO systems in a similar
way. It is also well-known that such representations are
much more convenient for a future controller design. As
far as the identification problem is concerned, a modified
predictor-based subspace identification (PBSID) (Chiuso,
2007) adapted with the propagator method (Munier and
Delisle, 1991; Merce`re and Lovera, 2007; Merce`re et al.,
2008; Merce`re and Bako, 2011), is used in order to estimate
directly a state-space realization of the system from the
measured input-output (I/O) data. Contrary to the classic
subspace algorithms (Van Overschee and De Moor, 1996;
Katayama, 2005; Verhaegen and Verdult, 2007), this tech-
nique does not give access to a fully-parametrized form but
leads to a state-space representation with a minimal num-
ber of parameters, even for MIMO systems. In addition,
thanks to the propagator method, the state vector can be
estimated directly from the available I/O data (Houtzager
et al., 2009). By doing so, a model linear with respect to
the parameters (LP) is obtained which makes the descrip-
tion of uncertainty areas easier. These domains are derived
from the analysis of particular quadratic criteria.
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
Copyright by the
International Federation of Automatic Control (IFAC)
4404
Page 2
The final objective is to get realistic uncertainty domains
that contain all kinds of stochastic disturbances. To reach
this goal, a hard unknown-but-bounded approach is con-
sidered. More particularly, an easy-tuning method is sug-
gested fixing the value of the required bound. Thus, no
specific assumption is made as far the noise is concerned.
In this paper, the method developed in (Farah et al.,
2010a) is applied to a novel wind turbine concept: the
smart rotor. This is a turbine equipped with additional
actuators and sensors to have a higher degree of freedom
for load mitigation (Van Wingerden et al., 2008, 2011;
Hulskamp et al., 2011). The controllers developed sofar
rely on estimated models without an uncertainty quantifi-
cation. The goal of this paper is to show the accuracy of
the obtain models.
The outline of this paper is as follows. In Section 2, the
main problem is stated. Section 3 is dedicated to the
estimation of the model parameters by using a particular
least-squares subspace algorithm. The uncertainty domain
determination problem is studied in Section 4. In Section 5,
the global technique performance is emphasized thanks to
experimental simulations. Section 6 concludes the paper.
2. PROBLEM FORMULATION AND NOTATIONS
Hereafter, we will consider a system with periodic distur-
bances, the shape of which is known a priori. It becomes
relevant if we are going to look at the actual dynamics of
the smart rotor. More precisely, the following state-space
model is used
x(t + 1) = Ax(t) +Bu(t) +Gd(t) + Fe(t) (1a)
y(t) = Cx(t) + e(t) (1b)
where x(t) ∈ Rnx , u(t) ∈ Rnu and y(t) ∈ Rny , are the
state, input and output vectors, respectively, e(t) ∈ Rny
denotes a zero mean white innovation process. Hereafter,
d(t) ∈ Rm contains the basis functions needed to express
the periodic disturbances. These periodic disturbances are
introduced because the system of interest, a two-bladed
smart rotor, significantly experiences periodic loads which
are multiples of the rotational speed. A, B, C, G and F
are the system, input, output, observer, periodic and noise
matrices respectively.
The main objective of this paper is to get reliable estimates
of the state-space matrices (A,B,G,C) combined with
accurate uncertainty parameter regions. As far as the
state-space matrices estimation is concerned, a new version
of the well-known PBSID algorithm is introduced (Chiuso,
2007). As shown hereafter, this algorithm gives access to
a consistent estimate of the state vector. This feature
will make easier the uncertainty domain description (see
Section 4).
The PBSID algorithm is based on the following predictor
form
x(t + 1) = A˜x(t) + B˜z(t) (2a)
y(t) = Cx(t) + e(t) (2b)
with A˜ = A − FC, B˜ = [B G F] and z⊤(t) =[
u⊤(t) d⊤(t) y⊤(t)
]⊤
. It corresponds to a straightforward
rewriting of the initial model (1). Hereafter, the following
notations will be used. More precisely, let us define the
stacked vector
up(t) =
[
u⊤(t− p) · · · u⊤(t− 1)
]⊤
where p is a user-defined integer such that p ≥ nx. In
a similar way, the stacked vectors yp(t) and zp(t) can
be obtained. The extended controllability matrix can be
defined as
K =
[
A˜p−1B˜ · · · A˜B˜ B˜
]
. (3)
We also define the observability matrix of the model (1),
the observability matrix of the model (2) and the Toeplitz
matrix
Γp =
[
(C)⊤ (CA)⊤ · · · (CAp−1)⊤
]⊤
Γ˜p =
[
(C)⊤ (CA˜)⊤ · · · (CA˜p−1)⊤
]⊤
H˜p =
0 0 · · · 0
CB˜ 0 · · · 0
...
...
. . .
...
CA˜p−2B˜ CA˜p−3B˜ · · · 0
.
With these definitions, the I/O behavior of the system (2)
is now given by (Verhaegen and Verdult, 2007)
yp(t) = Γ˜px(t− p) + H˜pzp(t) + ep(t). (4)
3. SUBSPACE-BASED LEAST-SQUARES
ALGORITHM (SBLS)
In the literature dedicated to subspace-based identification
(Van Overschee and De Moor, 1996; Katayama, 2005;
Verhaegen and Verdult, 2007), many algorithms try
(1) to get an accurate estimate xˆ of the state vector (via
suitable projections on available data sequences),
(2) to solve the following regression problem 1
min
A,B,C,D
∥∥∥∥
[
xˆ(t + 1)
y(t)
]
−
[
A B
C D
] [
x(t)
u(t)
]∥∥∥∥
2
F
(5)
in order to extract consistent estimates of the system
state-space matrices.
One of them is the well-known PBSID algorithm which is
based on a specific vector auto-regressive with exogenous
inputs model (Chiuso, 2007). This algorithm is interesting
because of its straightforward extension to the closed-loop
system identification problem. In this section, the basic
idea of the PBSID algorithm is used and adapted by
introducing a particular linear operator called the prop-
agator (Munier and Delisle, 1991; Merce`re and Lovera,
2007; Merce`re et al., 2008; Merce`re and Bako, 2011).
As a CCA-type algorithm (Bauer, 2005), the first step of
the developed algorithm herein consists in estimating the
state vector x. From Eq. (2), it is clear that (Peternell
et al., 1996; Knudsen, 2001; Chiuso, 2007)
x(t) = A˜px(t− p) +Kzp(t).
As with the PBSID algorithm, the key assumption of
our method is A˜j ≈ 0 for all j ≥ p. It can be shown
that, if the system (2) is uniformly exponential stable,
the approximation error can be made arbitrarily small
by making p large (Peternell et al., 1996; Knudsen, 2001;
1 ‖.‖F is the Frobenius norm.
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4405
that contain all kinds of stochastic disturbances. To reach
this goal, a hard unknown-but-bounded approach is con-
sidered. More particularly, an easy-tuning method is sug-
gested fixing the value of the required bound. Thus, no
specific assumption is made as far the noise is concerned.
In this paper, the method developed in (Farah et al.,
2010a) is applied to a novel wind turbine concept: the
smart rotor. This is a turbine equipped with additional
actuators and sensors to have a higher degree of freedom
for load mitigation (Van Wingerden et al., 2008, 2011;
Hulskamp et al., 2011). The controllers developed sofar
rely on estimated models without an uncertainty quantifi-
cation. The goal of this paper is to show the accuracy of
the obtain models.
The outline of this paper is as follows. In Section 2, the
main problem is stated. Section 3 is dedicated to the
estimation of the model parameters by using a particular
least-squares subspace algorithm. The uncertainty domain
determination problem is studied in Section 4. In Section 5,
the global technique performance is emphasized thanks to
experimental simulations. Section 6 concludes the paper.
2. PROBLEM FORMULATION AND NOTATIONS
Hereafter, we will consider a system with periodic distur-
bances, the shape of which is known a priori. It becomes
relevant if we are going to look at the actual dynamics of
the smart rotor. More precisely, the following state-space
model is used
x(t + 1) = Ax(t) +Bu(t) +Gd(t) + Fe(t) (1a)
y(t) = Cx(t) + e(t) (1b)
where x(t) ∈ Rnx , u(t) ∈ Rnu and y(t) ∈ Rny , are the
state, input and output vectors, respectively, e(t) ∈ Rny
denotes a zero mean white innovation process. Hereafter,
d(t) ∈ Rm contains the basis functions needed to express
the periodic disturbances. These periodic disturbances are
introduced because the system of interest, a two-bladed
smart rotor, significantly experiences periodic loads which
are multiples of the rotational speed. A, B, C, G and F
are the system, input, output, observer, periodic and noise
matrices respectively.
The main objective of this paper is to get reliable estimates
of the state-space matrices (A,B,G,C) combined with
accurate uncertainty parameter regions. As far as the
state-space matrices estimation is concerned, a new version
of the well-known PBSID algorithm is introduced (Chiuso,
2007). As shown hereafter, this algorithm gives access to
a consistent estimate of the state vector. This feature
will make easier the uncertainty domain description (see
Section 4).
The PBSID algorithm is based on the following predictor
form
x(t + 1) = A˜x(t) + B˜z(t) (2a)
y(t) = Cx(t) + e(t) (2b)
with A˜ = A − FC, B˜ = [B G F] and z⊤(t) =[
u⊤(t) d⊤(t) y⊤(t)
]⊤
. It corresponds to a straightforward
rewriting of the initial model (1). Hereafter, the following
notations will be used. More precisely, let us define the
stacked vector
up(t) =
[
u⊤(t− p) · · · u⊤(t− 1)
]⊤
where p is a user-defined integer such that p ≥ nx. In
a similar way, the stacked vectors yp(t) and zp(t) can
be obtained. The extended controllability matrix can be
defined as
K =
[
A˜p−1B˜ · · · A˜B˜ B˜
]
. (3)
We also define the observability matrix of the model (1),
the observability matrix of the model (2) and the Toeplitz
matrix
Γp =
[
(C)⊤ (CA)⊤ · · · (CAp−1)⊤
]⊤
Γ˜p =
[
(C)⊤ (CA˜)⊤ · · · (CA˜p−1)⊤
]⊤
H˜p =
0 0 · · · 0
CB˜ 0 · · · 0
...
...
. . .
...
CA˜p−2B˜ CA˜p−3B˜ · · · 0
.
With these definitions, the I/O behavior of the system (2)
is now given by (Verhaegen and Verdult, 2007)
yp(t) = Γ˜px(t− p) + H˜pzp(t) + ep(t). (4)
3. SUBSPACE-BASED LEAST-SQUARES
ALGORITHM (SBLS)
In the literature dedicated to subspace-based identification
(Van Overschee and De Moor, 1996; Katayama, 2005;
Verhaegen and Verdult, 2007), many algorithms try
(1) to get an accurate estimate xˆ of the state vector (via
suitable projections on available data sequences),
(2) to solve the following regression problem 1
min
A,B,C,D
∥∥∥∥
[
xˆ(t + 1)
y(t)
]
−
[
A B
C D
] [
x(t)
u(t)
]∥∥∥∥
2
F
(5)
in order to extract consistent estimates of the system
state-space matrices.
One of them is the well-known PBSID algorithm which is
based on a specific vector auto-regressive with exogenous
inputs model (Chiuso, 2007). This algorithm is interesting
because of its straightforward extension to the closed-loop
system identification problem. In this section, the basic
idea of the PBSID algorithm is used and adapted by
introducing a particular linear operator called the prop-
agator (Munier and Delisle, 1991; Merce`re and Lovera,
2007; Merce`re et al., 2008; Merce`re and Bako, 2011).
As a CCA-type algorithm (Bauer, 2005), the first step of
the developed algorithm herein consists in estimating the
state vector x. From Eq. (2), it is clear that (Peternell
et al., 1996; Knudsen, 2001; Chiuso, 2007)
x(t) = A˜px(t− p) +Kzp(t).
As with the PBSID algorithm, the key assumption of
our method is A˜j ≈ 0 for all j ≥ p. It can be shown
that, if the system (2) is uniformly exponential stable,
the approximation error can be made arbitrarily small
by making p large (Peternell et al., 1996; Knudsen, 2001;
1 ‖.‖F is the Frobenius norm.
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4405
Page 3
Chiuso, 2007). Under this assumption, the state x(t) is
approximately given by (Houtzager et al., 2009)
x(t) ≈ Kzp(t) (6)
where K is a matrix of unknown coefficients (see Eq.
(3)). To estimate this matrix, the following steps can be
considered. First of all,
y(t) ≈ CKzp(t) + e(t).
Hence, the matrices CK can be estimated by solving the
following linear problem
min
CK
‖y(t) −CKzp(t)‖2F .
For finite p, the solution of this linear problem will be
biased due to the approximation made in Eq. (6). In the
literature, several papers have studied the effect of the
window size and have proved the asymptotic properties of
the algorithms (if p → ∞ the bias disappears) (Knudsen,
2001; Chiuso, 2007). Then, with the approximation given
in Eq. (6), Eq. (4) can be written as
yp(t+ p) ≈ Γ˜pKzp(t) + H˜pzp(t+ p) + ep(t + p).
The product Γ˜pK is given by (Van Wingerden et al., 2008)
Γ˜pK =
CA˜p−1B˜ CA˜p−2B˜ · · · CB˜
CA˜pB˜ CA˜p−1B˜
. . . CA˜B˜
...
. . .
. . .
...
CA˜2p−2B˜ · · · · · · CA˜p−1B˜
.
Using the assumption A˜j ≈ 0 for all j ≥ p, this expression
can be approximated by the following upper block diagonal
matrix
Γ˜pK ≈
CA˜p−1B˜ CA˜p−2B˜ · · · CB˜
0 CA˜p−1B˜
. . . CA˜B˜
...
. . .
. . .
...
0 · · · · · · CA˜p−1B˜
. (7)
By introducing zeros in this matrix, the first block row
in Eq. (7) can be used to construct the other block rows.
Thus, knowing CK, the matrix Γ˜pK can be built.
In order to get estimates of (A,B,G,C) instead of
(A˜, B˜, G˜, C˜), the third step of the developed approach
consists in finding the relation between Γ˜pK and ΓpK.
As shown in (Houtzager et al., 2009), the following trans-
formation holds
ΓpK0
ΓpK1
...
ΓpKp−1
︸ ︷︷ ︸
ΓpK
=
Iny 0 · · · 0
CK Iny
. . .
...
...
. . .
. . . 0
CAp−2K CAp−1K · · · Iny
︸ ︷︷ ︸
T
Γ˜pK0
Γ˜pK1
...
Γ˜Kp−1
︸ ︷︷ ︸
Γ˜pK
which brings the observability matrix into the following
innovation form
ΓpK =
CA˜p−1B˜ CA˜p−2B˜ · · · CB˜
CAA˜p−1B˜ CAA˜p−2B˜ · · · CAB˜
...
. . .
. . .
...
CAp−1A˜p−1B˜ · · · · · · CAp−1B˜
.
The generators of the matrix T are not known, but, as
shown in (Dong et al., 2008), this transformation can be
performed as follows
ΓpKi = Γ˜pKi +
i−1∑
j=0
[
CA˜i−j−1K ΓpKi
]
with CA˜i−j−1K = CK(:, (p − i + j)(nu + ny) + nu +
(1 : ny)), i = 0 · · · p− 1 and where all the generators are
known.
The last step relies on a classic assumption in system iden-
tification and realization theory: the system is observable.
When this property is satisfied, because p ≥ nx, Γp has,
at least, nx linearly independent rows. Noting that Γp has
full column rank, a permutation matrix S ∈ Rnyp×nyp can
be introduced in order to reorder the rows of Γp and to
ensure that its first nx rows are linearly independent, i.e.,
SΓp =
[
Γ1
Γ2
]
where Γ1 is a matrix block containing a set of nx indepen-
dent rows and Γ2 is a matrix block containing the other
rows. Because Γ1 is a square full rank matrix, Γ2 can be
written as a linear combination of Γ1 or, similarly, there
is a unique operator P ∈ R(nyf−nx)×nx , named the prop-
agator (Munier and Delisle, 1991), such that Γ2 = PΓ1.
Then, it holds that
SΓpK =
[
Γ1
Γ2
]
K =
[
Inx
P
]
Γ1K
and
SΓpKzp(t) =
[
Inx
P
]
Γ1 Kzp(t)︸ ︷︷ ︸
≈x(t)
.
Consequently, knowing ΓpK and using Γ1 as a similarity
transformation, an estimate of the state vector in a new
state coordinate basis can be obtained
xˆ = Γ1Kzp(t)
from which an equivalent representation of the state-space
form (1) can be derived. Indeed, knowing y, u, d and xˆ,
the system matrices (A, B, G, C) can be estimated by
solving the following linear problem
min
A,B,G,C
∥∥∥∥∥
[
xˆ(t + 1)
y(t)
]
−
[
A B G
C 0 0
] [xˆ(t)
u(t)
d(t)
]∥∥∥∥∥
2
F
.
Because Γ1 is used as a similarity transformation, (A,
B, G, C) are realized in an observability canonical form
(Kailath, 1980; Merce`re and Bako, 2011). Thanks to this
property, the parameters of the estimated matrices are
invariant coefficients. This characteristic makes the de-
termination of the confidence intervals of the estimated
parameters much easier than with the classic subspace-
based algorithms.
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4406
approximately given by (Houtzager et al., 2009)
x(t) ≈ Kzp(t) (6)
where K is a matrix of unknown coefficients (see Eq.
(3)). To estimate this matrix, the following steps can be
considered. First of all,
y(t) ≈ CKzp(t) + e(t).
Hence, the matrices CK can be estimated by solving the
following linear problem
min
CK
‖y(t) −CKzp(t)‖2F .
For finite p, the solution of this linear problem will be
biased due to the approximation made in Eq. (6). In the
literature, several papers have studied the effect of the
window size and have proved the asymptotic properties of
the algorithms (if p → ∞ the bias disappears) (Knudsen,
2001; Chiuso, 2007). Then, with the approximation given
in Eq. (6), Eq. (4) can be written as
yp(t+ p) ≈ Γ˜pKzp(t) + H˜pzp(t+ p) + ep(t + p).
The product Γ˜pK is given by (Van Wingerden et al., 2008)
Γ˜pK =
CA˜p−1B˜ CA˜p−2B˜ · · · CB˜
CA˜pB˜ CA˜p−1B˜
. . . CA˜B˜
...
. . .
. . .
...
CA˜2p−2B˜ · · · · · · CA˜p−1B˜
.
Using the assumption A˜j ≈ 0 for all j ≥ p, this expression
can be approximated by the following upper block diagonal
matrix
Γ˜pK ≈
CA˜p−1B˜ CA˜p−2B˜ · · · CB˜
0 CA˜p−1B˜
. . . CA˜B˜
...
. . .
. . .
...
0 · · · · · · CA˜p−1B˜
. (7)
By introducing zeros in this matrix, the first block row
in Eq. (7) can be used to construct the other block rows.
Thus, knowing CK, the matrix Γ˜pK can be built.
In order to get estimates of (A,B,G,C) instead of
(A˜, B˜, G˜, C˜), the third step of the developed approach
consists in finding the relation between Γ˜pK and ΓpK.
As shown in (Houtzager et al., 2009), the following trans-
formation holds
ΓpK0
ΓpK1
...
ΓpKp−1
︸ ︷︷ ︸
ΓpK
=
Iny 0 · · · 0
CK Iny
. . .
...
...
. . .
. . . 0
CAp−2K CAp−1K · · · Iny
︸ ︷︷ ︸
T
Γ˜pK0
Γ˜pK1
...
Γ˜Kp−1
︸ ︷︷ ︸
Γ˜pK
which brings the observability matrix into the following
innovation form
ΓpK =
CA˜p−1B˜ CA˜p−2B˜ · · · CB˜
CAA˜p−1B˜ CAA˜p−2B˜ · · · CAB˜
...
. . .
. . .
...
CAp−1A˜p−1B˜ · · · · · · CAp−1B˜
.
The generators of the matrix T are not known, but, as
shown in (Dong et al., 2008), this transformation can be
performed as follows
ΓpKi = Γ˜pKi +
i−1∑
j=0
[
CA˜i−j−1K ΓpKi
]
with CA˜i−j−1K = CK(:, (p − i + j)(nu + ny) + nu +
(1 : ny)), i = 0 · · · p− 1 and where all the generators are
known.
The last step relies on a classic assumption in system iden-
tification and realization theory: the system is observable.
When this property is satisfied, because p ≥ nx, Γp has,
at least, nx linearly independent rows. Noting that Γp has
full column rank, a permutation matrix S ∈ Rnyp×nyp can
be introduced in order to reorder the rows of Γp and to
ensure that its first nx rows are linearly independent, i.e.,
SΓp =
[
Γ1
Γ2
]
where Γ1 is a matrix block containing a set of nx indepen-
dent rows and Γ2 is a matrix block containing the other
rows. Because Γ1 is a square full rank matrix, Γ2 can be
written as a linear combination of Γ1 or, similarly, there
is a unique operator P ∈ R(nyf−nx)×nx , named the prop-
agator (Munier and Delisle, 1991), such that Γ2 = PΓ1.
Then, it holds that
SΓpK =
[
Γ1
Γ2
]
K =
[
Inx
P
]
Γ1K
and
SΓpKzp(t) =
[
Inx
P
]
Γ1 Kzp(t)︸ ︷︷ ︸
≈x(t)
.
Consequently, knowing ΓpK and using Γ1 as a similarity
transformation, an estimate of the state vector in a new
state coordinate basis can be obtained
xˆ = Γ1Kzp(t)
from which an equivalent representation of the state-space
form (1) can be derived. Indeed, knowing y, u, d and xˆ,
the system matrices (A, B, G, C) can be estimated by
solving the following linear problem
min
A,B,G,C
∥∥∥∥∥
[
xˆ(t + 1)
y(t)
]
−
[
A B G
C 0 0
] [xˆ(t)
u(t)
d(t)
]∥∥∥∥∥
2
F
.
Because Γ1 is used as a similarity transformation, (A,
B, G, C) are realized in an observability canonical form
(Kailath, 1980; Merce`re and Bako, 2011). Thanks to this
property, the parameters of the estimated matrices are
invariant coefficients. This characteristic makes the de-
termination of the confidence intervals of the estimated
parameters much easier than with the classic subspace-
based algorithms.
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4406
Page 4
4. UNCERTAINTY DOMAIN DESCRIPTION
As soon as the user has access to the estimated state-
space matrices, the uncertainty domain determination
problem can be considered. More precisely, the uncertainty
domain description of the coefficients ai, b
j
i , g
k
i and c
l
j ,
i ∈ [0, nx − 1], j ∈ [1, nu], k ∈ [2, ny] and l ∈ [1, ny] (i.e.
the state-space matrices A, B, G and C parameters) can
be performed. The basic idea of the approach introduced
hereafter consists in studying particular ellipsoidal iso-
level curves by adapting a bounded error method for which
only a bounded residual hypothesis is required (Baron
et al., 2001). This original method has been developed
by assuming that the model of the system is described
via a linear regression. In order to use efficiently this
technique in the framework considered in this paper, the
model (1) must be transformed into a linear regression.
More precisely, knowing y, u, d and xˆ, it is straightforward
that
vec
[
x(t + 1)⊤ y(t)⊤
]
︸ ︷︷ ︸
s(t)
=
[
Inx+ny ⊗
[
x(t)⊤ u(t)⊤ d(t)⊤
]⊤]
︸ ︷︷ ︸
Φ⊤(t)
vec
A⊤ C⊤
B⊤ 0⊤
G⊤ 0⊤
︸ ︷︷ ︸
θ
+ vec
([
Fe(t)
e(t)
]⊤)
︸ ︷︷ ︸
v(t)
. (8)
In (Ljung, 1999a), a technique is developed to describe
the uncertainty domain by using the statistical properties
of the prediction error (PE) method. In this technique,
the input sequence is assumed to be persistently exciting
(Ljung, 1999a) and the system to be stable. Using these
assumptions, the PE identification procedure delivers an
unbiased identified parameter vector. More precisely, this
vector is asymptotically a random vector with gaussian
distribution, a mean equal to the true parameter vector θ0
and covariance equal to λ( 1N
∑N
t=1 Φ
⊤(t)Φ(t))−1 where λ
is the noise variance. In a different way,
θ̂ ∈ N
(
θ0, λ( 1
N
N∑
t=1
Φ⊤(t)Φ(t))−1
)
where θ̂ is an arbitrary estimate of θ. Consequently, the
uncertainty domain for the true parameter vector θ0 can
be given by
(θ̂ − θ0)⊤( 1
N
N∑
t=1
Φ⊤(t)Φ(t))(θ̂ − θ0) < αλ (9)
with α is a user-defined probability level. This is an
ellipsoidal domain centered at θ0 whose shape is given
by 1N
∑N
t=1 Φ
⊤(t)Φ(t). This approach assumes that θ̂ is
asymptotically a random vector with gaussian distribution
and mean value θ0. Furthermore, the description of this
confidence domain depends on the prior knowledge of the
noise variance.
Because our objective is to characterize the uncertainty
domain without any assumption on the noise but only the
information contained in the residuals, the idea presented
in (Ljung, 1999a) is adapted hereafter. To do so, let
us go back to Eq. (8). Then, by analyzing this linear
regression, it is obvious that the least-squares estimation
of this parameter vector can be obtained by minimizing
the following cost function (Ljung, 1999a, Appendix II)
J(θ) = 1
N
N∑
t=1
ε⊤(t, θ)ε(t, θ)
where ε(t, θ) = s(t) −Φ⊤(t)θ. Furthermore, it is easy to
see that (Ljung, 1999a)
J(θ) = 1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂)
︸ ︷︷ ︸
computable residuals
+(θ̂ − θ)⊤RΦ(θ̂ − θ)︸ ︷︷ ︸
ellipsoid centered on θ̂
where RΦ = 1N
∑N
t=1Φ
⊤(t)Φ(t). This specific decompo-
sition shows that the cost function J(θ) is made up of
two terms: the mean square of the residuals (which is
computable a posteriori) as well as an ellipsoid centered
on θˆ. Noticing this decomposition, the uncertainty domain
determination problem can be solved by fixing “efficiently”
the level J(θ) so that the true parameter vector θ0 is
ensured to be included in the confidence ellipsoid. In order
to find this consistent level, let us observe that, when
θ = θ0, we have
(θ̂ − θ0)⊤RΦ(θ̂ − θ0) = J(θ0)−
1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂)).
Thus, by choosing a threshold M as follows
M = J(θ0)− 1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂),
it is guaranteed that the real parameter vector is located
on the border of the uncertainty ellipsoidal domain. By
construction, this bound M is a function of the unknown
parameter vector θ0. Thus, in practice, we can only choose
a bound M˜
M˜ > J(θ0)− 1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂).
As suggested in (Ljung, 1999a), this bound M˜ could be
chosen from an estimate of the noise signal v because, by
construction, it can be easily seen that
M˜ = J(θ0) = 1
N
N∑
t=1
v⊤(t)v(t) (10)
which leads to the following inequality
J(θ0)− 1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂) < 1
N
N∑
t=1
v⊤(t)v(t).
This way of fixing the bound M˜ requires reliable knowl-
edge of the noise v. In order to circumvent this difficulty,
it is suggested fixing the size of the ellipsoid domain by us-
ing exclusively the information contained in the residuals
ε(t, θ̂). In order to see the link between the residuals and
the noise, let us decompose ε(t, θ̂) as follows
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4407
As soon as the user has access to the estimated state-
space matrices, the uncertainty domain determination
problem can be considered. More precisely, the uncertainty
domain description of the coefficients ai, b
j
i , g
k
i and c
l
j ,
i ∈ [0, nx − 1], j ∈ [1, nu], k ∈ [2, ny] and l ∈ [1, ny] (i.e.
the state-space matrices A, B, G and C parameters) can
be performed. The basic idea of the approach introduced
hereafter consists in studying particular ellipsoidal iso-
level curves by adapting a bounded error method for which
only a bounded residual hypothesis is required (Baron
et al., 2001). This original method has been developed
by assuming that the model of the system is described
via a linear regression. In order to use efficiently this
technique in the framework considered in this paper, the
model (1) must be transformed into a linear regression.
More precisely, knowing y, u, d and xˆ, it is straightforward
that
vec
[
x(t + 1)⊤ y(t)⊤
]
︸ ︷︷ ︸
s(t)
=
[
Inx+ny ⊗
[
x(t)⊤ u(t)⊤ d(t)⊤
]⊤]
︸ ︷︷ ︸
Φ⊤(t)
vec
A⊤ C⊤
B⊤ 0⊤
G⊤ 0⊤
︸ ︷︷ ︸
θ
+ vec
([
Fe(t)
e(t)
]⊤)
︸ ︷︷ ︸
v(t)
. (8)
In (Ljung, 1999a), a technique is developed to describe
the uncertainty domain by using the statistical properties
of the prediction error (PE) method. In this technique,
the input sequence is assumed to be persistently exciting
(Ljung, 1999a) and the system to be stable. Using these
assumptions, the PE identification procedure delivers an
unbiased identified parameter vector. More precisely, this
vector is asymptotically a random vector with gaussian
distribution, a mean equal to the true parameter vector θ0
and covariance equal to λ( 1N
∑N
t=1 Φ
⊤(t)Φ(t))−1 where λ
is the noise variance. In a different way,
θ̂ ∈ N
(
θ0, λ( 1
N
N∑
t=1
Φ⊤(t)Φ(t))−1
)
where θ̂ is an arbitrary estimate of θ. Consequently, the
uncertainty domain for the true parameter vector θ0 can
be given by
(θ̂ − θ0)⊤( 1
N
N∑
t=1
Φ⊤(t)Φ(t))(θ̂ − θ0) < αλ (9)
with α is a user-defined probability level. This is an
ellipsoidal domain centered at θ0 whose shape is given
by 1N
∑N
t=1 Φ
⊤(t)Φ(t). This approach assumes that θ̂ is
asymptotically a random vector with gaussian distribution
and mean value θ0. Furthermore, the description of this
confidence domain depends on the prior knowledge of the
noise variance.
Because our objective is to characterize the uncertainty
domain without any assumption on the noise but only the
information contained in the residuals, the idea presented
in (Ljung, 1999a) is adapted hereafter. To do so, let
us go back to Eq. (8). Then, by analyzing this linear
regression, it is obvious that the least-squares estimation
of this parameter vector can be obtained by minimizing
the following cost function (Ljung, 1999a, Appendix II)
J(θ) = 1
N
N∑
t=1
ε⊤(t, θ)ε(t, θ)
where ε(t, θ) = s(t) −Φ⊤(t)θ. Furthermore, it is easy to
see that (Ljung, 1999a)
J(θ) = 1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂)
︸ ︷︷ ︸
computable residuals
+(θ̂ − θ)⊤RΦ(θ̂ − θ)︸ ︷︷ ︸
ellipsoid centered on θ̂
where RΦ = 1N
∑N
t=1Φ
⊤(t)Φ(t). This specific decompo-
sition shows that the cost function J(θ) is made up of
two terms: the mean square of the residuals (which is
computable a posteriori) as well as an ellipsoid centered
on θˆ. Noticing this decomposition, the uncertainty domain
determination problem can be solved by fixing “efficiently”
the level J(θ) so that the true parameter vector θ0 is
ensured to be included in the confidence ellipsoid. In order
to find this consistent level, let us observe that, when
θ = θ0, we have
(θ̂ − θ0)⊤RΦ(θ̂ − θ0) = J(θ0)−
1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂)).
Thus, by choosing a threshold M as follows
M = J(θ0)− 1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂),
it is guaranteed that the real parameter vector is located
on the border of the uncertainty ellipsoidal domain. By
construction, this bound M is a function of the unknown
parameter vector θ0. Thus, in practice, we can only choose
a bound M˜
M˜ > J(θ0)− 1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂).
As suggested in (Ljung, 1999a), this bound M˜ could be
chosen from an estimate of the noise signal v because, by
construction, it can be easily seen that
M˜ = J(θ0) = 1
N
N∑
t=1
v⊤(t)v(t) (10)
which leads to the following inequality
J(θ0)− 1
N
N∑
t=1
ε⊤(t, θ̂)ε(t, θ̂) < 1
N
N∑
t=1
v⊤(t)v(t).
This way of fixing the bound M˜ requires reliable knowl-
edge of the noise v. In order to circumvent this difficulty,
it is suggested fixing the size of the ellipsoid domain by us-
ing exclusively the information contained in the residuals
ε(t, θ̂). In order to see the link between the residuals and
the noise, let us decompose ε(t, θ̂) as follows
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4407
Page 5
ε(t, θ̂) = s(t)−Φ⊤(t)θ̂ = v(t) −Φ⊤(t)(θ0 − θ̂).
Thus, the residuals contain information related to the
noise as well as the modelling error θ0 − θ̂. Furthermore,
by assuming that the regressor Φ and the noise v are
uncorrelated, it is straightforward that 2
E
{
ε⊤(t, θ̂)ε(t, θ̂)
}
= E
{
v⊤(t)v(t) + ∆θ⊤Φ(t)Φ⊤(t)∆θ
}
(11)
with ∆θ = θ0 − θ̂. This equality has two interesting
properties. First, when θ̂ is an unbiased estimate of θ0,
the last term of Eq. (11) disappears and we have a reliable
estimate of the mean square of the noise. Second, when
the estimation is biased, the mean square of the residuals
contains the error due to the noise v as well as the
modeling error ∆θ. In order to take into account these
two interesting properties, the bound M˜ is chosen equal
to
M˜ =
1
N
max
N∑
t=1
(
ε⊤(t, θ̂)ε(t, θ̂)
)
. (12)
5. EXPERIMENTAL EXAMPLE
In this section, we present a prototyped smart rotor
used to show the feasibility of the ‘smart’ rotor concept
with the advanced control (for more details see (van
Wingerden, 2008)). The ‘smart’ rotor that we use for our
experimental validation is a rotating, two-bladed rotor
equipped with trailing-edge flaps. Each blade is equipped
with two distinct trailing-edge flaps in order to enable its
use for future research. The main focus for the controller
is to suppress the bending modes of the two individual
blades and consequently the two actuators on one blade are
used together as one actuator by applying the same control
signal. For the same reasons, two strain sensors are applied
in the root located on the central axis and at the leading
edge of the blade, respectively. The experimental setup
mainly consists of the following components: wind tunnel,
blade, actuators & sensors, and real-time environment (for
more details see (Van Wingerden et al., 2011)). Because
the final goal for this rotating smart rotor is to reduce the
fatigue loads, two macro fiber composite (MFC) patches
are adhered to the root to measure the high strains
associated with the first flap wise bending mode. So, it is
interesting to identify the model between the trailing-edge
flaps on the two different blades u, and the corresponding
MFC sensors y. The identification procedure is carried out
with a rotation speed equal to 430 RPM. The order of the
system is estimated to be equal to 12. To assess the quality
of the models, the Best Fit (% BFT) quality of the model
(Ljung, 1999a) and the Percent Variance Accounted For
(% VAF) (Verhaegen and Verdult, 2007) are used. This
criterions are defined as
BFT = 100%×max
(
1− ‖y − ŷ‖2‖y − ym‖2
, 0
)
(13)
V AF = 100%×max
(
1− var(y − ŷ)
var(y)
, 0
)
(14)
2 E{•} stands for the mathematical expectation of the signal •.
Blade 1 Blade 2
Algorithm % BFT % VAF % BFT % VAF
SBLS 73.93 93.20 72.15 96.01
PO-MOESP 66.57 88.82 71.72 95.43
Table 1. % BFT and % VAF for the different
identification methods.
where ŷ is the simulated output and ym is the mean
value of y. In Table 1, the algorithm introduced in this
paper (named SBLS) is compared with the PO-MOESP
algorithm (Verhaegen and Verdult, 2007) by looking at the
% VAF and % BFT. Table 1 shows that the SBLS method
gives better results than the PO-MOESP method.
Because the system contains a large number of parameters,
it is difficult to display all the estimated parameters as
well as the corresponding uncertainty domains. It is more
convenient to present the frequency response uncertainty
domain by using the following equation
cov(∆G∆G∗) =
δĜ
δθ R
δĜ
δθ
∗
(15)
with G(ω) = C(jωI −A)−1B. By using Ljung’s method
(Ljung, 1997), R = λ(Φ⊤Φ)−1 and by using the uncer-
tainty method presented in Section 4, R = M˜(Φ⊤Φ)−1.
The confidence bound on the frequency response is pre-
sented in Fig. 1. These curves show that the developed
method is more conservative than Ljung’s one. Neverthe-
less, this method takes into account the modeling error.
These uncertainty regions present a great interest in a
robust control objective.
10−2 100
−100
−90
−80
−70
−60
−50
−40
input1
o
u
tp
ut
1
10−2 100
−80
−70
−60
−50
input2
10−2 100
−70
−65
−60
−55
−50
−45
input1
o
u
tp
ut
2
10−2 100
−60
−55
−50
−45
input2
Fig. 1. estimated (blue solid line) frequency responses with
confidence intervals given by the developed method
(red dashed lines) and Ljung’s 3σ confidence method
(black dotted line).
6. CONCLUSION
In this paper, the problem of the state-space coefficients
uncertainty domain determination for MIMO LTI state-
space systems has been tackled. More precisely, a par-
ticular least-squares subspace-based method has been in-
troduced in order to identify the state-space model in
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4408
Thus, the residuals contain information related to the
noise as well as the modelling error θ0 − θ̂. Furthermore,
by assuming that the regressor Φ and the noise v are
uncorrelated, it is straightforward that 2
E
{
ε⊤(t, θ̂)ε(t, θ̂)
}
= E
{
v⊤(t)v(t) + ∆θ⊤Φ(t)Φ⊤(t)∆θ
}
(11)
with ∆θ = θ0 − θ̂. This equality has two interesting
properties. First, when θ̂ is an unbiased estimate of θ0,
the last term of Eq. (11) disappears and we have a reliable
estimate of the mean square of the noise. Second, when
the estimation is biased, the mean square of the residuals
contains the error due to the noise v as well as the
modeling error ∆θ. In order to take into account these
two interesting properties, the bound M˜ is chosen equal
to
M˜ =
1
N
max
N∑
t=1
(
ε⊤(t, θ̂)ε(t, θ̂)
)
. (12)
5. EXPERIMENTAL EXAMPLE
In this section, we present a prototyped smart rotor
used to show the feasibility of the ‘smart’ rotor concept
with the advanced control (for more details see (van
Wingerden, 2008)). The ‘smart’ rotor that we use for our
experimental validation is a rotating, two-bladed rotor
equipped with trailing-edge flaps. Each blade is equipped
with two distinct trailing-edge flaps in order to enable its
use for future research. The main focus for the controller
is to suppress the bending modes of the two individual
blades and consequently the two actuators on one blade are
used together as one actuator by applying the same control
signal. For the same reasons, two strain sensors are applied
in the root located on the central axis and at the leading
edge of the blade, respectively. The experimental setup
mainly consists of the following components: wind tunnel,
blade, actuators & sensors, and real-time environment (for
more details see (Van Wingerden et al., 2011)). Because
the final goal for this rotating smart rotor is to reduce the
fatigue loads, two macro fiber composite (MFC) patches
are adhered to the root to measure the high strains
associated with the first flap wise bending mode. So, it is
interesting to identify the model between the trailing-edge
flaps on the two different blades u, and the corresponding
MFC sensors y. The identification procedure is carried out
with a rotation speed equal to 430 RPM. The order of the
system is estimated to be equal to 12. To assess the quality
of the models, the Best Fit (% BFT) quality of the model
(Ljung, 1999a) and the Percent Variance Accounted For
(% VAF) (Verhaegen and Verdult, 2007) are used. This
criterions are defined as
BFT = 100%×max
(
1− ‖y − ŷ‖2‖y − ym‖2
, 0
)
(13)
V AF = 100%×max
(
1− var(y − ŷ)
var(y)
, 0
)
(14)
2 E{•} stands for the mathematical expectation of the signal •.
Blade 1 Blade 2
Algorithm % BFT % VAF % BFT % VAF
SBLS 73.93 93.20 72.15 96.01
PO-MOESP 66.57 88.82 71.72 95.43
Table 1. % BFT and % VAF for the different
identification methods.
where ŷ is the simulated output and ym is the mean
value of y. In Table 1, the algorithm introduced in this
paper (named SBLS) is compared with the PO-MOESP
algorithm (Verhaegen and Verdult, 2007) by looking at the
% VAF and % BFT. Table 1 shows that the SBLS method
gives better results than the PO-MOESP method.
Because the system contains a large number of parameters,
it is difficult to display all the estimated parameters as
well as the corresponding uncertainty domains. It is more
convenient to present the frequency response uncertainty
domain by using the following equation
cov(∆G∆G∗) =
δĜ
δθ R
δĜ
δθ
∗
(15)
with G(ω) = C(jωI −A)−1B. By using Ljung’s method
(Ljung, 1997), R = λ(Φ⊤Φ)−1 and by using the uncer-
tainty method presented in Section 4, R = M˜(Φ⊤Φ)−1.
The confidence bound on the frequency response is pre-
sented in Fig. 1. These curves show that the developed
method is more conservative than Ljung’s one. Neverthe-
less, this method takes into account the modeling error.
These uncertainty regions present a great interest in a
robust control objective.
10−2 100
−100
−90
−80
−70
−60
−50
−40
input1
o
u
tp
ut
1
10−2 100
−80
−70
−60
−50
input2
10−2 100
−70
−65
−60
−55
−50
−45
input1
o
u
tp
ut
2
10−2 100
−60
−55
−50
−45
input2
Fig. 1. estimated (blue solid line) frequency responses with
confidence intervals given by the developed method
(red dashed lines) and Ljung’s 3σ confidence method
(black dotted line).
6. CONCLUSION
In this paper, the problem of the state-space coefficients
uncertainty domain determination for MIMO LTI state-
space systems has been tackled. More precisely, a par-
ticular least-squares subspace-based method has been in-
troduced in order to identify the state-space model in
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4408
Page 6
a fixed coordinate basis. This identification method is
made up of two steps. The first leads to a consistent
estimate of the state vector thanks to the propagator
method in a user-defined coordinate basis. Then, the state
matrices (A,B,G,C) are estimated using a least-squares
algorithm. The uncertainty domains are described using
a bounded error approach. This approach assumes that
the residual is bounded and no stochastic assumption on
the noise is made. This method has been finally applied
in order to identify the model of a smart rotor. These
experimental results have highlighted the performance of
the developed method.
REFERENCES
C. Baron, T. Poinot, and J. Trigeassou. Determination
of parametric uncertainty domains using least-squares
technique and bounded errors. Methods and Models in
Automation and Control, 2:1003–1008, 2001.
D. Bauer. Asymptotic properties of subspace estimators.
Automatica, 41:359–376, 2005.
A. Chiuso. The role of vector autoregressive modeling in
predictor-based subspace identification. Automatica, 43:
1034–1048, 2007.
D. de Vries. Identification of model uncertainty for control
design. PhD thesis, Delft University of Technology,
Delft, The Netherlands, 1994.
J. Dong, M. Verhaegen, and E. Holweg. Closed-loop sub-
space predictive control for fault tolerant MPC design.
In Proceedings of the 17th World Congress IFAC, Seoul,
Korea, February 2008.
W. Farah, G. Merce`re, and T. Poinot. System identifica-
tion and uncertainty domain determination: a subspace-
based approach. In Proceedings of the American Control
Conference, Baltimore, Maryland, USA, June 2010a.
W. Farah, G. Merce`re, J. W. van Wingerden, and
T. Poinot. Quantification of model uncertainty for a
state-space system. In Proceedings of the International
Symposium on Mathematical Theory of Networks and
Systems, Budapest, Hungary, July 2010b.
L. Giarre´ and M. Milanese. Model quality evaluation
in H2 identification. IEEE Transactions on Automatic
Control, 42:691–698, 1997.
L. Giarre´, M. Milanese, and M. Taragna. H∞ identification
and model quality evaluation. IEEE Transactions on
Automatic Control, 42:188–199, 1997.
G. Goodwin, M. Gevers, and B. Ninness. Quantifying the
error in estimated transfer functions with application to
model error selection. IEEE Transactions on Automatic
Control, 37:913–928, 1992.
I. Houtzager, J. W. van Wingerden, and M. Verhaegen.
Fast-array recursive closed-loop subspace model identi-
fication. In Proceedings of the 15th IFAC Symposium on
System Identification, Saint-Malo, France, July 2009.
A. Hulskamp, J. W. Van Wingerden, T. Barlas, H. Cham-
pliaud, G. Van Kuik, H. Bersee, and M. Verhaegen.
Dynamic load mitigation experiments on a scale model
of a wind turbine with a smart rotor. To appear in Wind
Energy, 2011.
T. Kailath. Linear Systems. Prentice Hall, Engelwood
Cliffs, 1980.
T. Katayama. Subspace methods for system identification.
Springer Verlag, 2005.
T. Knudsen. Consistency analysis of subspace identifi-
cation methods based on a linear regression approach.
Automatica, 37:81–89, 2001.
L. Ljung. Identification model validation and control. In
Proceedings of the IEEE Conference on Decision and
Control, San Diego, California, USA, December 1997.
L. Ljung. System identification. Theory for the user.
Prentice Hall, Upper Saddle River, 2nd edition, 1999a.
L. Ljung. Model validation and model error modeling.
In Proceedings of the Astro¨m Symposium on Control,
Lund, Sweden, August 1999b.
G. Merce`re and L. Bako. Parameterization and identifi-
cation of multivariable state-space systems: a canonical
approach. Accepted provisionally as Regular Paper in
Automatica, 2011.
G. Merce`re and M. Lovera. Convergence analysis of
instrumental variable recursive subspace identification
algorithms. Automatica, 43:1377–1386, 2007.
G. Merce`re, L. Bako, and S. Lecœuche. Propagator-based
methods for recursive subspace model identification.
Signal Processing, 88:468–491, 2008.
J. Munier and G. Delisle. Spatial analysis using new prop-
erties of the cross spectral matrix. IEEE Transactions
on Signal Processing, 39:746–749, 1991.
B. Ninness and G. Goodwin. Estimation of model quality.
Automatica, 31:1771–1797, 1995.
K. Peternell, W. Sherrer, and M. Deistler. Statistical
analysis of novel subspace identification methods. Signal
Processing, 52:161–177, 1996.
W. Reinelt, G. Garulli, and L. Ljung. Comparing different
approaches to model error modeling in robust identifi-
cation. Automatica, 38:787–803, 2001.
P. Van Overschee and B. De Moor. Subspace identification
for linear systems. theory, implementation, applications.
Kluwer Academic Publishers, 1996.
J. W. van Wingerden. Control of wind turbines with smart
rotors: proof of concept and LPV subspace identification.
PhD thesis, Delft University of Technology, Delft, The
Netherlands, 2008.
J. W. Van Wingerden, A. Hulskamp, T. Barlas, G. Van
Kuik, D. Molenaar, H. Bersee, and M. Verhaegen. On
the proof of concept of a smart wind turbine rotor blade
for load alleviation. Wind Energy, 11:265–280, 2008.
J. W. Van Wingerden, A. Hulskamp, T. Barlas, I. Houtza-
ger, H. Bersee, G. Van Kuik, and M. Verhaegen. Two-
degree-of-freedom active vibration control of a proto-
typed smart rotor. IEEE Transactions on Control Sys-
tems Technology, 19:284–296, 2011.
M. Verhaegen and V. Verdult. Filtering and system
identification: a least squares approach. Cambridge
University Press, 2007.
E. Walter, J. Norton, H. Piet-Lahanier, and M. Milanese.
Bounding Approaches to System Identification. Perseus
Publishing, 1996.
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4409
made up of two steps. The first leads to a consistent
estimate of the state vector thanks to the propagator
method in a user-defined coordinate basis. Then, the state
matrices (A,B,G,C) are estimated using a least-squares
algorithm. The uncertainty domains are described using
a bounded error approach. This approach assumes that
the residual is bounded and no stochastic assumption on
the noise is made. This method has been finally applied
in order to identify the model of a smart rotor. These
experimental results have highlighted the performance of
the developed method.
REFERENCES
C. Baron, T. Poinot, and J. Trigeassou. Determination
of parametric uncertainty domains using least-squares
technique and bounded errors. Methods and Models in
Automation and Control, 2:1003–1008, 2001.
D. Bauer. Asymptotic properties of subspace estimators.
Automatica, 41:359–376, 2005.
A. Chiuso. The role of vector autoregressive modeling in
predictor-based subspace identification. Automatica, 43:
1034–1048, 2007.
D. de Vries. Identification of model uncertainty for control
design. PhD thesis, Delft University of Technology,
Delft, The Netherlands, 1994.
J. Dong, M. Verhaegen, and E. Holweg. Closed-loop sub-
space predictive control for fault tolerant MPC design.
In Proceedings of the 17th World Congress IFAC, Seoul,
Korea, February 2008.
W. Farah, G. Merce`re, and T. Poinot. System identifica-
tion and uncertainty domain determination: a subspace-
based approach. In Proceedings of the American Control
Conference, Baltimore, Maryland, USA, June 2010a.
W. Farah, G. Merce`re, J. W. van Wingerden, and
T. Poinot. Quantification of model uncertainty for a
state-space system. In Proceedings of the International
Symposium on Mathematical Theory of Networks and
Systems, Budapest, Hungary, July 2010b.
L. Giarre´ and M. Milanese. Model quality evaluation
in H2 identification. IEEE Transactions on Automatic
Control, 42:691–698, 1997.
L. Giarre´, M. Milanese, and M. Taragna. H∞ identification
and model quality evaluation. IEEE Transactions on
Automatic Control, 42:188–199, 1997.
G. Goodwin, M. Gevers, and B. Ninness. Quantifying the
error in estimated transfer functions with application to
model error selection. IEEE Transactions on Automatic
Control, 37:913–928, 1992.
I. Houtzager, J. W. van Wingerden, and M. Verhaegen.
Fast-array recursive closed-loop subspace model identi-
fication. In Proceedings of the 15th IFAC Symposium on
System Identification, Saint-Malo, France, July 2009.
A. Hulskamp, J. W. Van Wingerden, T. Barlas, H. Cham-
pliaud, G. Van Kuik, H. Bersee, and M. Verhaegen.
Dynamic load mitigation experiments on a scale model
of a wind turbine with a smart rotor. To appear in Wind
Energy, 2011.
T. Kailath. Linear Systems. Prentice Hall, Engelwood
Cliffs, 1980.
T. Katayama. Subspace methods for system identification.
Springer Verlag, 2005.
T. Knudsen. Consistency analysis of subspace identifi-
cation methods based on a linear regression approach.
Automatica, 37:81–89, 2001.
L. Ljung. Identification model validation and control. In
Proceedings of the IEEE Conference on Decision and
Control, San Diego, California, USA, December 1997.
L. Ljung. System identification. Theory for the user.
Prentice Hall, Upper Saddle River, 2nd edition, 1999a.
L. Ljung. Model validation and model error modeling.
In Proceedings of the Astro¨m Symposium on Control,
Lund, Sweden, August 1999b.
G. Merce`re and L. Bako. Parameterization and identifi-
cation of multivariable state-space systems: a canonical
approach. Accepted provisionally as Regular Paper in
Automatica, 2011.
G. Merce`re and M. Lovera. Convergence analysis of
instrumental variable recursive subspace identification
algorithms. Automatica, 43:1377–1386, 2007.
G. Merce`re, L. Bako, and S. Lecœuche. Propagator-based
methods for recursive subspace model identification.
Signal Processing, 88:468–491, 2008.
J. Munier and G. Delisle. Spatial analysis using new prop-
erties of the cross spectral matrix. IEEE Transactions
on Signal Processing, 39:746–749, 1991.
B. Ninness and G. Goodwin. Estimation of model quality.
Automatica, 31:1771–1797, 1995.
K. Peternell, W. Sherrer, and M. Deistler. Statistical
analysis of novel subspace identification methods. Signal
Processing, 52:161–177, 1996.
W. Reinelt, G. Garulli, and L. Ljung. Comparing different
approaches to model error modeling in robust identifi-
cation. Automatica, 38:787–803, 2001.
P. Van Overschee and B. De Moor. Subspace identification
for linear systems. theory, implementation, applications.
Kluwer Academic Publishers, 1996.
J. W. van Wingerden. Control of wind turbines with smart
rotors: proof of concept and LPV subspace identification.
PhD thesis, Delft University of Technology, Delft, The
Netherlands, 2008.
J. W. Van Wingerden, A. Hulskamp, T. Barlas, G. Van
Kuik, D. Molenaar, H. Bersee, and M. Verhaegen. On
the proof of concept of a smart wind turbine rotor blade
for load alleviation. Wind Energy, 11:265–280, 2008.
J. W. Van Wingerden, A. Hulskamp, T. Barlas, I. Houtza-
ger, H. Bersee, G. Van Kuik, and M. Verhaegen. Two-
degree-of-freedom active vibration control of a proto-
typed smart rotor. IEEE Transactions on Control Sys-
tems Technology, 19:284–296, 2011.
M. Verhaegen and V. Verdult. Filtering and system
identification: a least squares approach. Cambridge
University Press, 2007.
E. Walter, J. Norton, H. Piet-Lahanier, and M. Milanese.
Bounding Approaches to System Identification. Perseus
Publishing, 1996.
Preprints of the 18th IFAC World Congress
Milano (Italy) August 28 - September 2, 2011
4409
Sign up today - FREE
Mendeley saves you time finding and organizing research. Learn more
- All your research in one place
- Add and import papers easily
- Access it anywhere, anytime
Start using Mendeley in seconds!
Readership Statistics
5 Readers on Mendeley
by Discipline
20% Engineering
by Academic Status
60% Ph.D. Student
20% Student (Bachelor)
20% Associate Professor
by Country
40% France
20% Brazil
20% Hong Kong


