Quantification of model uncertainty for a state-space system
Wafa Farah, Guillaume Mercère, Thierry Poinot and Jan-Willem van Wingerden
Abstract— In this communication, the uncertainty domain
determination problem for multi-input multi-output systems de-
scribed with a linear time-invariant state-space representation
is adressed. The developed method is based on a two-step ap-
proach. The first step consists in estimating the nominal model
using a particular least-squares subspace algorithm. Then, the
uncertainty domains are described by using a bounded error
approach. Simulations are used to highlight the performance
of the method.
I. INTRODUCTION
One of the main objectives of control-oriented identifica-
tion is to estimate models that are suitable for robust control
design techniques [1]. For this purpose, system identification
must give not only a nominal model, but also a reliable
estimate of the uncertainty associated with the model. To
reach this goal, two steps are generally necessary. The first
one deals with the model parameters estimation. This model
is required to understand, to control or to improve the
system functioning. The second step consists in designing the
control law from the model parameters. However, because
the model is only a system approximation, 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 estimated model parameters
uncertainties. Thus, the estimated model uncertainties must
be well-described.
In system identification theory, the uncertainty domain
description is mainly based on prior assumptions about
noise and unmodeled dynamics [2], [1]. The first works
assume that the disturbances acting on the system are random
variables realizations [3], [4]. Because the information on
measurement noises are not often available and difficult to
verify [5], a different characterization way can be considered
[6]. This approach is mainly based on deterministic hypothe-
ses, i.e. on the assumption that the residuals are unknown-
but-bounded [7], [8], [9]. This basic idea has given rise to
a number of techniques usually addressed as bounded error
or set membership identification [10], [11], [12], [13]. 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
W. Farah, G. Mercère and T. Poinot are with the Uni-
versity of Poitiers, Laboratoire d’Automatique et Informatique
Industrielle, 40, avenue du recteur Pineau, 86022 Poitiers,
France wafa.farah@etu.univ-poitiers.fr,
guillaume.mercere@univ-poitiers.fr,
thierry.poinot@univ-poitiers.fr
Jan-Willem van Wingerden is with the University of Delft, Delft
Center of System and Control, Delft, 2628 CD The Netherlands.
J.W.vanWingerden@tudelft.nl
and the noise affecting the data) which makes the bound
determination quite difficult. To tackle this problem, some
works [4], [14], [15] propose to estimate the error modeling.
In [16], the bounded error is determined through the analysis
of particular iso-level curves which leads to less conservative
bounds. In other words, this bounding approach uses only an
assumption on the residual bound and not a stochastic noise
assumption. The identification and uncertainty domain de-
scription methods developed in [16] are restricted to single-
input single-output (SISO) system represented by transfer
function. So, it will be interesting to extend this method to
multi-input multi-output (MIMO) system.
Because it is more convenient to use the state-space
representation in the MIMO case, a subspace method will
be used to estimate directly a state-space realization of the
system from the measured input-output (I/O). The originality
of the approach proposed hereafter, is to determine the state-
space parameters uncertainty using bounded error approach.
To reach this goal, a particular subspace-based method,
named the propagator method [17], [18], will be used to
estimate directly a state-space realization of the system from
the measured I/O data. Contrary to the classic subspace
algorithms [19], [20], [21], this technique does not give
access to a fully-parametrized form but leads to a state-space
representation with a minimal number of parameters, even
for MIMO systems. This method gives also a state-space
model estimate in a user-defined state-space coordinates
basis. In addition, thanks to the propagator method, the
state vector can be estimated [22]. So, a model linear with
respect to the parameters (LP) is obtained which makes the
description of uncertainty areas easier. These domains are
derived from the analysis of particular quadratic criteria in
the optimum vicinity. The final objective is to get realistic
uncertainty domains that contain all kinds of stochastic
disturbances. Thus, a hard unknown-but-bounded approach
is considered to reach this goal. More particularly, an easy-
tuning method is proposed to fix the value of the required
bound and no assumption is made on the noise.
The outline of this paper is as follows. In Section II
the main problem is stated. Section III is dedicated to the
system parameters estimation using a particular least-squares
subspace algorithm. The uncertainty domain determination
problem is studied in Section IV. In Section V, the global
technique performance is emphasized thanks to numerical
simulations. Section VI concludes the paper.
Proceedings of the 19th International Symposium on Mathematical Theory of Networks and Systems – MTNS 2010 • 5–9 July, 2010 • Budapest, Hungary
ISBN 978-963-311-370-7 1007

II. PROBLEM FORMULATION
Consider the following linear time-invariant state-space
system
x(t+ 1) = Ax(t) +Bu(t) (1a)
y(t) = Cx(t) +Du(t) + v(t). (1b)
where u(t) ∈ Rnu , y(t) ∈ Rny , x(t) ∈ Rnx and v(t) ∈
Rny are respectively the input, the output, the state and the
output-noise vectors. (A,B,C,D) are the system matrices
relatively to a certain state-space coordinate basis. The model
(1) is the so-called output-error model [23]. The system order
is assumed to be known a priori1.
Notice that no particular assumption is made concerning
the properties of v. However, it is assumed hereafter that
the model residuals, i.e., the difference between the sim-
ulated and measured outputs, are bounded. This particular
framework can be justified by noticing that, in many prac-
tical cases, residuals are related to complex dynamics that
cannot be captured by the plant model. When the residuals
contain some deterministic parts, it turns out that the classic
stationary stochastic process assumption is no more suitable.
On the contrary, deterministic constraint such as “bounded
in magnitude residuals” may be more efficient [10].
The considered identification problem can be stated as
follows: given realizations {u(t)}Nt=1 and {y(t)}
N
t=1 of the
input and output processes generated by a system of the
form (1) on a finite but sufficiently wide time horizon N ,
the goal of the developed method is to estimate the matrices
(A,B,C,D) represented in a particular state-space form and
characterize the state-space parameters precision. To reach
this goal, two linked problems are solved. Firstly, the system
matrices (1) are identified using a particular least-squares
subspace algorithm (see Section III). Then, the estimated
parameters uncertainty domains are characterized through the
analysis of particular iso-level curves and an unknown-but-
bounded approach (see Section IV).
III. SUBSPACE-BASED LEAST-SQUARES
ALGORITHM
The proposed uncertainty domain method assumes that the
model is LP. For this reason, a particular subspace-based
least-squares algorithm is adopted to identify the model
parameters. This algorithm has the particularity to estimate
the state vector which allows to get an LP model. To reach
this goal, some definitions are introduced in §III-A. Then,
the subspace-based least-squares algorithm is introduced in
§III-B.
1This assumption is not really strong because many algorithms, mainly
based on SVD and information criteria [20], are now available to get a
reliable estimate of this parameter. Thus, for LTI systems, the system order
can be estimated beforehand.
A. Data equation
Let us define
up(t) =
[
u⊤(t− p) · · · u⊤(t− 1)
]⊤
K =
[
Ap−1B · · · AB B
]
Γp =
[
(C)⊤ (CA)⊤ · · · (CAp−1)⊤
]⊤
Hp =


D 0 · · · 0
CB D · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
CAp−2B CAp−3B · · · D


where p is a user-defined integer such that p ≥ nx. With
these definitions, the I/O behavior of the system (1) is now
given by [21]
yp(t) = Γpx(t− p) +Hpup(t) + vp(t). (2)
B. Identification algorithm
Because the uncertainty domain determination problem
is easier when the model is LP, the state vector is firstly
estimated. This vector is identified thanks to the particular
property of the propagator method. Then the (A,B,C,D)
matrices are estimated using a least-squares algorithm.
The starting point of this method is the following well-
known relation [24], [25], [26]
x(t) = Apx(t− p) +Kup(t).
The key assumption in this method is that we assume that
Aj ≈ 0 for all j ≥ p. It can be shown that if the system (1)
is uniformly exponential stable, the approximation error can
be made arbitrarily small by making p large [24][27][25].
Under this assumption, the state x(t) is approximately
given by [22]
x(t) ≈ Kup(t) (3)
and the output behavior can be written as
y(t) ≈ CKup(t) +Du(t) + v(t).
Hence, the matrices CK can be estimated by solving the
following linear problem
min
CK
∥∥∥∥y(t)−
[
CK D
] [up(t)
u(t)
]∥∥∥∥
2
F
.
For finite p, the solution of this linear problem will be biased
due to the approximation made in (3). 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) [24][27][25].
With the approximation given in (3), we can rewrite (2)
as
yp(t+ p) ≈ ΓpKup(t) +Hpup(t+ p) + vp(t+ p).
The product ΓpK is given by [28]
ΓpK =


CAp−1B CAp−2B · · · CB
CApB CAp−1B
.
.
. CAB
.
.
.
.
.
.
.
.
.
.
.
.
CA2p−2B · · · · · · CAp−1B


.
W. Farah et al. • Quantification of Model Uncertainty for a State-Space System
1008

Using the assumption Aj ≈ 0 for all j ≥ p, this expression
can be approximated by the following upper block diagonal
matrix
ΓpK ≈


CAp−1B CAp−2B · · · CB
0 CAp−1B
.
.
. CAB
.
.
.
.
.
.
.
.
.
.
.
.
0 · · · · · · CAp−1B


. (4)
By introducing zeros in this matrix, the first block row
in (4) can be used to construct the other block rows. So,
knowing CK, the matrix ΓpK can be constructed.
Because our objective is to get an estimate of x(t), the
propagator method will be used [18][17]. Firstly, thanks to
the estimate of ΓpK, the signal
q = ΓpKup(t) ≈ Kx(t)
can be simulated. Then, assuming that the system (1) is
observable, Γp has got, at least, nx linearly independent
rows. Introducing a permutation matrix S ∈ Rnyp×nyp such
that the extended observability matrix can be reordered as
follows
q =
[
q1
q2
]
= SΓpKup =
[
Γ1
Γ2
]
Kup =
[
Inx
P
]
Γ1Kup
where Γ1 is a matrix block containing a set of nx indepen-
dent rows, Γ2 is a matrix block containing the other rows
and P is a unique operator named the propagator [17].
It holds that
xˆ(t) ≈ Γ1Kup(t). (5)
So, The state vector can be estimated in a user-defined
coordinate basis.
Knowing y, u and xˆ, the system matrices (A, B, C, D)
can be estimated by solving the following linear problem
min
ABCD
∥∥∥∥
[
x
y
]
−
[
A B
C D
] [
x
u
]∥∥∥∥
2
F
.
Using Γ1 as a similarity transformation, the state-space
matrices satisfy
A =
[
Inx
P
]
(2 : nx + 1, :) (6a)
=


0 1 0 · · · 0
0 0 1 · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 · · · 1
−a0 −a1 −a2 · · · −anx−1


(6b)
c1 =
[
1 0 · · · 0
]
cj =
[
Inx
P
]
((j − 1)i+ 1, :) for j ∈ [2, ny] (6c)
with aj , j ∈ [0, nx − 1], the coefficients of the characteristic
polynomial of A [29]. Knowing A and c, the matrices B
and D can be estimated by using a least squares regression
[21]
B̂, D̂ = argminB,D
∑
t [y(t)
−
[
u⊤(t)⊗ Iny
∑t−1
k=0 u⊤(k)⊗ ĉÂ
t−k−1
] [
vec(D)
vec(B)
]]
where vec(.) stands for the vectorization operator and ⊗ is
the Kronecker product [30].
Thanks to the propagator method, (A,B,C,D) are ex-
pressed in LP form. So, the proposed method to describe the
uncertainty domain can be used to determin the (A,B,C,D)
parameters uncertainty.
IV. UNCERTAINTY DOMAINS DESCRIPTION
As soon as the user has access to the estimated state-space
matrices, the uncertainty domains determination problem can
be considered. More precisely, the uncertainty domain de-
scription of the coefficients ai, bji , cki and dlj , i ∈ [0, nx − 1],
j ∈ [1, nu], k ∈ [2, ny] and l ∈ [1, ny] (the state-space matri-
ces A, B, C and D parameters) is performed. The developed
approach basic idea is to study particular ellipsoidal iso-level
curves by adapting a bounded error method in which only
bounded residual hypothesis is required [16]. This original
method has been developed assuming that the system model
is described via a linear regression.
For this reason, knowing y, u and xˆ, the model (1) can
be transformed as follow
vec
[
x(t+ 1)⊤ y(t)⊤
]
︸ ︷︷ ︸
yM
=
[
Inx+ny ⊗
[
x(t)⊤ u(t)⊤
]⊤]
︸ ︷︷ ︸
Φ⊤(t)
vec
[
A⊤ C⊤
B⊤ D⊤
]
︸ ︷︷ ︸
θ
+vec(v(t)).
Remark 1 The first nx−1 rows of A and the first row of C
are made up of 0 and 1. For this reason, θ can be restricted
as
θ˜ = vec
[
A(nx, :)⊤ C(j, :)⊤
B⊤ D⊤
]
for j ∈ [2, ny] .
yM and Φ can be deduced in a compatible way with θ˜.
It is obvious that the least-squares estimates of these
matrices can be obtained by minimizing the cost function
[1, Appendix II]2
J(θ̂) = 12
(
yM −ΨM θ̂
)⊤ (
yM −ΨM θ̂
)
= 12ε
⊤ε
2In this method, We don’t take into account the uncertainty of the state
sequence estimate.
Proceedings of the 19th International Symposium on Mathematical Theory of Networks and Systems – MTNS 2010 • 5–9 July, 2010 • Budapest, Hungary
1009

where θ̂ is an arbitrary estimate of θ,
yM =


vec[x(2)⊤ y(1)⊤]
.
.
.
vec[x(M + 1)⊤ y(M)]


and ΨM =


Φ⊤(1)
.
.
.
Φ⊤(M)

 .
Then, assuming that RΨ = Ψ⊤MΨM is invertible, it is easy
to see that
J(θ̂) = 12y
⊤
M
(
I−ΨM
(
Ψ⊤MΨM
)−1Ψ⊤M
)
yM (7)
+ 12
(
θ̂ − θls
)⊤
Ψ⊤MΨM
(
θ̂ − θls
)
= Jmin +
1
2dθ̂
⊤RΨdθ̂
with θls =
(
Ψ⊤MΨM
)−1Ψ⊤MyM and dθ̂ the variation of
the estimate θ̂ around θls. This relation shows that J(θ̂)
has a unique minimum at the least-squares solution θls and
the first term of the lhs of Eq. (7) is the minimum value.
Furthermore, if we introduce V (θ̂) = J(θ̂) − Jmin, it is
obvious that
V (θ̂) = 12dθ̂
⊤RΨdθ̂ (8)
is an ellipsoid centered in θls whose main directions are
given by RΨ.
The goal of the developed method is to determine an
uncertainty domain D such that θ ∈ D. Here, D will be
an ellipsoidal surface depending on a user-defined criterion
level JD. This threshold must be chosen such that the system
parameters θ is surely included in D without leading to a
too conservative solution. This level is obtained hereafter by
using an analogy with the stochastic framework. To reach
this goal, this link will be presented firstly in gaussian case
and then in general case.
Under the Gaussian assumptions (no modeling error and
zero-mean white Gaussian output noise), it is well-known
that [1]
1
2 (θls − θ)
⊤Ψ⊤MΨM (θls − θ) = n2σ2
where σ2 is the noise variance and n ∈ R+. This is again an
ellipsoid with the same shape as the one given by (8). This
analogy leads to choose the level JD as follows
JD − Jmin = n2σ2 = ℓ2.
In the general case, when the disturbance v is not Gaussian
but only zero-mean and modeling error is considered, we can
say that
−nσ < v < nσ with n ≥ 3 and ℓ2 = n2σ2. (9)
So, an ellipsoid domain D including all the possible values
of θ is defined by
JD = Jmin + ℓ2 with ℓ = nσ. (10)
Unfortunately, this relation depends on the prior knowl-
edge of σ2. In order to circumvent this difficulty, a bounded
error approach is used3. To estimate this bound, symbolized
hereafter by ℓ, the first idea could be to measure the upper
bound of the noise acting on the system when the system
is not excited. Practically, this methodology cannot involve
the modeling error brought out when the system is excited.
Then, it is proposed to use the information contained in the
residuals. Indeed, the following result can be proved
Jmin =
1
2v
⊤
(
I−ΨM
(
Ψ⊤MΨM
)−1Ψ⊤M
)
v = 12ε
⊤ε
where v is the system disturbance and ε = yM − yˆM (θ) are
the residuals. Thus, a good noise effect approximation can
be deduced from ε and we fix ℓ as follows
ℓ = max {|ε(t)|} .
ℓ depend only on the residual and not on the noise acting on
the output.
V. SIMULATION EXAMPLE
In order to show the performance of the method described
beforehand, the following state-space matrices are used
A =


0.2 0 −1
1 0.3 5
−2 −0.4 −0.6

 B =


1 0
0 2
1 −1


C =
[
5 0 1
−3 1 1
]
.
These matrices can be rewritten as
A =


0 1 0
0 0 1
0.1 0.2 −0.1

 B =


7.3 0.2
−0.3 −6.9
−14.6 −2.0


C =
[
1.0 0 0
1.0 −0.3 0.6
]
by using
Γ1 =


5 0 1
−1 −0.4 −5.6
10.6 2.12 2.36


as similarity transformation.
The input signal is a pseudo random binary sequences
(PRBS) of length 1000. A Monte Carlo simulation of size
1000 is carried out. The output noise signal v is a zero-
mean white Gaussian noise such that the signal to noise ratio
(SNR) equals 20dB. To satisfy the assumption Aj ≈ 0 for all
j ≥ p, the past index is choosen as p = 20. The least-squares
method is applied to estimate the system matrices. The
results of this identification method are presented in Table
I. This table shows that the subspace-based least-squares
algorithm gives a good (A,B,C) matrices estimation.
3In this method, the only assumption is to have bounded residuals.
W. Farah et al. • Quantification of Model Uncertainty for a State-Space System
1010

True parameters Estimated parameters
−a0 0.1640 0.1640 ± 2.60 10−7
−a1 0.2400 0.2400 ± 3 10−7
−a2 −0.1000 −0.1000 ± 3.30 10−7
c20 1.0406 1.0406 ± 3.2 10−7
c21 −0.3068 −0.3068 ± 4.3 10−7
c22 0.6707 0.6706 ± 6.4 10−7
b10 7.3777 7.3775 ± 8.57 10−5
b20 0.2006 0.2007 ± 8.93 10−5
b11 −0.3933 −0.3936 ± 9.77 10−5
b21 −6.9678 −6.9682 ± 8.53 10−5
b12 −14.6091 −14.6093 ± 9.02 10−5
b22 −2.0070 −2.0070 ± 8.09 10−5
TABLE I
ESTIMATED PARAMETERS USING THE LEAST-SQUARES METHOD.
Identification results are plotted also in Figure 1. Firstly,
these three graphs show that the estimates obtained using
the least-squares method are accurate. Indeed, in each plot,
the black cross (+), symbolizing the mean value of the
estimated parameters, almost matches the red one (×) which
corresponds to the real parameters.
Figure 1 shows the system parameters (×), the mean value
of the estimated parameters (+) and the estimated parameters
(∗) calculated during the 1000 realizations of the Monte
Carlo simulation for (−a0,−a1), (−a2, b10) and (−a2, b10).
To assess the quality of these uncertainty domains and
to get rid of the drawing of all the uncertainty domains, a
failure rate measure is used. This failure rate is defined as the
percentage of realizations for which the system parameters
are outside of the ellipsoid D centered in the estimated
parameters vector, (see black disc (•) in Fig. 1). The failure
ratio equals 0.8 % for the couple (−a0,−a1), 2 % for
(−a2, b10) and 6.7 % for (−a2, b10). These values show that
the way the threshold is fixed leads to reliable uncertainty
domains.
VI. CONCLUSIONS
In this paper, the problem of the state-space coefficients
uncertainty domain determination for MIMO LTI state-space
systems is solved. Because the ordinary subspace method
provide an estimation of the system matrices up to a similar-
ity transformation, a particular least-squares method is used.
This identification method contains two steps. The first leads
to estimate the state vector thanks to the propagator method
in a user-define coordinate basis. Then, (A,B,C,D) are
estimated using a least-squares algorithm. The uncertainty
domains are described using a bounded error approach.
This approach supposes that the residual is bounded and no
stochastic assumption on the noise is made. The experimental
results have emphasized the reliability of the developed
method.
REFERENCES
[1] L. Ljung, System identification: theory for the user (2nd edition).
Prentice Hall, 1999.
0.162 0.163 0.164 0.165
0.239
0.24
0.241
0.242
−a0
−
a 1
0.162 0.163 0.164 0.165 0.166
7.35
7.36
7.37
7.38
7.39
7.4
7.41
−a0
b1 0
1.039 1.04 1.041 1.042
−0.308
−0.307
−0.306
−0.305
c20
c2 1
Fig. 1. Level surfaces of the cost function J(θ¯). The real parameters are
symbolized by a red cross (×), the estimated parameters by a blue cross
(∗) and the mean value of the estimated parameters by a black cross (+).
Black discs (•) are finally the failure draws.
[2] G. Goodwin, M. Gevers, and B. Ninness, “Quantifying the error in
estimated transfer functions with application to model order selection,”
IEEE Transactions on Automatic Control, vol. 37, pp. 913–928, 1992.
[3] L. Ljung, “Identification model validation and control,” in Proceed-
ings of the 36th Conference on Decision and Control, San Diego,
California, USA, December 1997.
[4] ——, “Model validation and model error modeling,” in Proceedings
of the Aström Symposium on Control, Lund, Sweden, August 1999.
[5] D. de Vries, “Identification of model uncertainty for control design,”
Ph.D. dissertation, Delft University of Technology, Delft, The Nether-
lands, 1994.
[6] B. Ninness and G. Goodwin, “Estimation of model quality,” Automat-
ica, vol. 31, pp. 1771–1797, 1995.
[7] L. Giarré and M. Milanese, “Model quality evaluation in H2 identifica-
tion,” IEEE Transactions on Automatic Control, vol. 42, pp. 691–698,
1997.
[8] L. Giarré, M. Milanese, and M. Taragna, “H∞ identification and
model quality evaluation,” IEEE Transactions on Automatic Control,
vol. 42, pp. 188–199, 1997.
Proceedings of the 19th International Symposium on Mathematical Theory of Networks and Systems – MTNS 2010 • 5–9 July, 2010 • Budapest, Hungary
1011

[9] W. Farah, G. Mercère, and T. Poinot, “System identification and
uncertainty domain determination: a subspace-based approach,” in
American Control Conference, Baltimore, USA, July 2010.
[10] E. Walter, J. Norton, H. Piet−Lahanier, and M. Milanese, Bounding
approaches to system identification. Plenum Press, 1996.
[11] G. Belforte, B. Dona, and V. Cerone, “Parameter estimation algorithms
for a set−membership description of uncertainty,” Automatica, vol. 26,
pp. 887–898, 1990.
[12] Y. Huang and S. Gollamudi, “Set-membership identification for adap-
tive equalization,” in Proceedings of the 38th Midwest Symposium on
Circuits and Systems, Rio de Janeiro, Brazil, August 1995.
[13] N. Ramdani and P. Poignet, “Experimental parallel robot dynamic
model evaluation with set membership estimation,” in Proceedings
of the 14th IFAC Symposium on System Identification, Newcastle,
Australia, March 2006.
[14] C. Baron, “Contribution à l’estimation de paramètres physiques à
l’aide du modèles d’ordre réduit,” Ph.D. dissertation, Université de
Poitiers, France, 2003.
[15] W. Reinelt, G. Garulli, and L. Ljung, “Comparing different approaches
to model error modeling in robust identification,” Automatica, vol. 38,
pp. 787–803, 2001.
[16] C. Baron, T. Poinot, and J.-C. Trigeassou, “Determination of paramet-
ric uncertainty domains using least-squares technique and bounded
errors,” Methods and Models in Automation and Control, vol. 2, pp.
1003–1008, August 2001.
[17] J. Munier and G. Y. Delisle, “Spatial analysis using new properties of
the cross spectral matrix,” IEEE Transactions on Signal Processing,
vol. 39, pp. 746–749, 1991.
[18] G. Mercère, L. Bako, and S. Lecœuche, “Propagator-based methods for
recursive subspace model identification,” Signal Processing, vol. 88,
pp. 468–491, 2008.
[19] P. V. Overschee and B. D. Moor, Subspace identification for linear
systems theory, implementation, applications. Kluwer Academic
Publishers, 1996.
[20] T. Katayama, Subspace methods for system identification. Springer-
Verlag, 2005.
[21] M. Verhaegen and V. Verdult, Filtering and system identification: a
least-squares approach. Cambridge University Press, 2007.
[22] I. Houtzager, J. W. van Wingerden, and M. Verhaegen, “Fast-array
recursive closed-loop subspace model identification,” in Proceedings
of the 15th IFAC symposium identification, Saint-Malo, France, July
2009.
[23] M. Verhaegen and P. Dewilde, “Subspace model identification part 2:
analysis of the elementry output error state space model identification
algorithm,” International Journal of Control, vol. 56, pp. 1211–1241,
1992.
[24] T. Knudsen, “Consistency analysis of subspace identification methods
based on a linear regression approach,” Automatica, vol. 37, pp. 81–89,
2001.
[25] A. Chiuso, “The role of vector auto regressive modeling in predictor
based subspace identification,” Automatica, vol. 43, pp. 1034–1048,
2007.
[26] K. Peternell, W. Sherrer, and M. Deistler, “Statistical analysis of novel
subspace identification methods,” Signal Processing, vol. 52, pp. 161–
177, 1996.
[27] A. Chiuso and G. Picci, “Consistency analysis of some closed-loop
subspace identification methods,” Automatica, vol. 41, pp. 371–391,
2005.
[28] J. W. van Wingerden, “Control of wind turbines with smart rptors:
proof of concept and lpv subspace identification,” Ph.D. dissertation,
Delft University, Delft, Netherlands, 2008.
[29] L. Bako, G. Mercère, and S. Lecœuche, “On-line structured subspace
identification with application to switched linear systems,” Interna-
tional Journal of Control, vol. 82, pp. 1496–1515, 2009.
[30] J. W. Brewer, “Kronecker product and matrix calculus in system
theory,” IEEE Transactions on Circuits and Systems, vol. 25, pp. 772–
781, 1978.
W. Farah et al. • Quantification of Model Uncertainty for a State-Space System
1012