System identification and uncertainty domain determination: a subspace-based approach
Available from
Guillaume Mercère's profile on Mendeley.
Page 1
System identification and uncertainty domain determination: a subspace-based approach
System identification and uncertainty domain determination: a
subspace-based approach
Wafa Farah, Guillaume Mercère and Thierry Poinot
Abstract—In this paper, the problem of uncertainty domain
determination for multi-input multi-output systems described
with linear time-invariant state-space representation is tackled.
The developed method combines a subspace-based identification
algorithm and a bounded-error approach. To estimate the
system state-space matrices, the propagator method is applied.
This method is characterized by the use of least-squares
minimization algorithms and the fact that the estimated state-
space realization is described using a fixed structure with a
minimal number of parameters. The problem of uncertainty
domain description is then solved from the analysis of particular
iso-level curves in the vicinity of the optimum. A hard error
bounding approach is finally considered to choose a threshold
required by the developed method.
I. INTRODUCTION
Generally speaking, in the automatic control framework,
the goal of system identification is to provide a system
model accurate for control design [1]. To reach this goal,
two steps are generally necessary. The first one deals with
the estimation of the model parameters. 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 uncertainties of the estimated model
parameters. 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 ran-
dom variables realizations [3], [4]. Mainly based on time-
domain representations [5], these first developments have
quickly shown their limits, particularly for robust control
[6]. Furthermore, in some practical cases, the assumption
of random white noise is too conservative or difficult to
verify a priori [7]. For this reason, other approaches have
been developed. Several of them rely on deterministic hy-
potheses, i.e. on non-probabilistic or hard error bounding
approaches [8], [9]. In these techniques, the identification
error is unknown but bounded. This basic idea has given
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
rise to a number of techniques usually addressed as bounded-
error or set membership identification [10]. Although many
set membership methods have been successfully applied on
real systems [11], [12], the main drawback of this approach
is its dependence on the way the bound is determined a
priori. Notice indeed that the error comes from two sources
(the unmodeled dynamics and the noise affecting the data)
which makes the bound determination quite difficult. To deal
with this problem, some works [4] propose to estimate the
error modeling. In this case, the bounded-error is determined
through the analysis of particular iso-level curves which
leads to less conservative bounds. The identification and
uncertainty domain description methods developed in [13]
are restricted to single-input single-output (SISO) system
represented by transfer function. The main goal of this
communication is to extend the basic idea of this approach
to multi-input multi-output (MIMO) system. In order to deal
in a similar way with SISO and MIMO systems and because
it is much more convenient for control design, state-space
representations will be used. As far as the identification
problem is concerned, a particular subspace-based method,
named the propagator method [14], [15], will be used to
estimate directly a state-space realization of the system
from the measured input-output (I/O) data. Contrary to the
classic subspace algorithms [16], [17], [18], this technique
does not give access to fully-parametrized form but leads
to a state-space representation with a minimal number of
parameters, even for MIMO systems. It is quite obvious that
the problem of uncertainty domain determination is easier
when the number of estimated parameters is constant and
minimal. Thanks to the fixed and minimal state-space ma-
trices structure provided by the identification algorithm, the
model is written into an input/output (I/O) form linear in the
parameters (LP) 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
error bounding approach is considered to reach this goal.
More particularly, an easy-tuning method is proposed to fix
the value of the required bound.
The outline of this paper is as follows. In Section II the
notations, the main problem and the general assumptions
are stated. Section III is dedicated to the system parameters
estimation using the propagator method. 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
2010 American Control Conference
Marriott Waterfront, Baltimore, MD, USA
June 30-July 02, 2010
FrC07.3
978-1-4244-7427-1/10/$26.00 ©2010 AACC 6501
subspace-based approach
Wafa Farah, Guillaume Mercère and Thierry Poinot
Abstract—In this paper, the problem of uncertainty domain
determination for multi-input multi-output systems described
with linear time-invariant state-space representation is tackled.
The developed method combines a subspace-based identification
algorithm and a bounded-error approach. To estimate the
system state-space matrices, the propagator method is applied.
This method is characterized by the use of least-squares
minimization algorithms and the fact that the estimated state-
space realization is described using a fixed structure with a
minimal number of parameters. The problem of uncertainty
domain description is then solved from the analysis of particular
iso-level curves in the vicinity of the optimum. A hard error
bounding approach is finally considered to choose a threshold
required by the developed method.
I. INTRODUCTION
Generally speaking, in the automatic control framework,
the goal of system identification is to provide a system
model accurate for control design [1]. To reach this goal,
two steps are generally necessary. The first one deals with
the estimation of the model parameters. 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 uncertainties of the estimated model
parameters. 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 ran-
dom variables realizations [3], [4]. Mainly based on time-
domain representations [5], these first developments have
quickly shown their limits, particularly for robust control
[6]. Furthermore, in some practical cases, the assumption
of random white noise is too conservative or difficult to
verify a priori [7]. For this reason, other approaches have
been developed. Several of them rely on deterministic hy-
potheses, i.e. on non-probabilistic or hard error bounding
approaches [8], [9]. In these techniques, the identification
error is unknown but bounded. This basic idea has given
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
rise to a number of techniques usually addressed as bounded-
error or set membership identification [10]. Although many
set membership methods have been successfully applied on
real systems [11], [12], the main drawback of this approach
is its dependence on the way the bound is determined a
priori. Notice indeed that the error comes from two sources
(the unmodeled dynamics and the noise affecting the data)
which makes the bound determination quite difficult. To deal
with this problem, some works [4] propose to estimate the
error modeling. In this case, the bounded-error is determined
through the analysis of particular iso-level curves which
leads to less conservative bounds. The identification and
uncertainty domain description methods developed in [13]
are restricted to single-input single-output (SISO) system
represented by transfer function. The main goal of this
communication is to extend the basic idea of this approach
to multi-input multi-output (MIMO) system. In order to deal
in a similar way with SISO and MIMO systems and because
it is much more convenient for control design, state-space
representations will be used. As far as the identification
problem is concerned, a particular subspace-based method,
named the propagator method [14], [15], will be used to
estimate directly a state-space realization of the system
from the measured input-output (I/O) data. Contrary to the
classic subspace algorithms [16], [17], [18], this technique
does not give access to fully-parametrized form but leads
to a state-space representation with a minimal number of
parameters, even for MIMO systems. It is quite obvious that
the problem of uncertainty domain determination is easier
when the number of estimated parameters is constant and
minimal. Thanks to the fixed and minimal state-space ma-
trices structure provided by the identification algorithm, the
model is written into an input/output (I/O) form linear in the
parameters (LP) 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
error bounding approach is considered to reach this goal.
More particularly, an easy-tuning method is proposed to fix
the value of the required bound.
The outline of this paper is as follows. In Section II the
notations, the main problem and the general assumptions
are stated. Section III is dedicated to the system parameters
estimation using the propagator method. 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
2010 American Control Conference
Marriott Waterfront, Baltimore, MD, USA
June 30-July 02, 2010
FrC07.3
978-1-4244-7427-1/10/$26.00 ©2010 AACC 6501
Page 2
paper.
II. PROBLEM FORMULATION
Consider the following LTI state-space model
x(t + 1) = Ax(t) +Bu(t) (1a)
y(t) = Cx(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) are the system matrices relatively
to a certain coordinate state-space basis. It is assumed that
{x(t)}, {u(t)}, {y(t)} and {v(t)}, t ∈ Z, are all ergodic
and (weakly) stationary stochastic processes [1]. The system
order is assumed1 to be known a priori.
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 , estimate
the matrices (A,B,C) represented in a particular state-space
form and characterize the estimated parameters precision.
We start by making the following assumptions: i) the
input process {u(t)} is an ergodic and (weakly) stationary
process that is persistently exciting, ii) the output noise
{v(t)} is a zero mean noise with a finite variance statistically
uncorrelated with the input {u(t)}, iii) the system is stable
and minimal. To reach the goal mentioned above, two linked
problems are solved. Firstly, the system matrices (1) are
identified using a particular subspace-based algorithm (see
Section III). Then, the uncertainty domains of the estimated
parameters are characterized through the analysis of partic-
ular iso-level curves and an error bounded approach (see
Section IV).
III. THE PROPAGATOR METHOD
A. Basic idea of the method
Subspace-based model identification (SMI) is a particular
field of experimental modeling which has reached an inter-
esting maturity from now on [19], [16], [20], [21], [17], [18].
SMI methods are attractive because state-space realizations
can be estimated directly from I/O data without using non-
linear optimization but by applying robust linear algebra
tools such as the RQ factorization and the singular value
decomposition (SVD) [22]. Most of the SMI algorithms
lead to minimal fully-parametrized state-space realizations
[16], [17], [18]. Interesting from a numerical point of view,
the use of fully-parametrized model structure clearly over-
parametrizes the model and results in identifiability loss [23,
Chapter 2]. Furthermore, by applying a SVD, the person
using the classic subspace algorithms has difficulties to fix
the basis in which the state-space models are estimated. In
fact, this choice is realized during the reduction step (when
the system order is estimated) but the reproducibility of
the basis is not ensured when such a tool is employed. In
1This assumption is not really strong because many algorithms, mainly
based on SVD and information criteria [17], are now available to get a
reliable estimate of this parameter. Thus, for LTI systems, the system order
can be estimated beforehand.
this paper, a solution is proposed to give access to models
described in a particular structured form where the number
of elements in the state-space matrices is nx(nu + ny), i.e.
the minimal number of parameters of any MIMO system
parametrization [24].
The key problem of this identification method is the
consistent estimation, from measured I/O samples, of the
extended observability matrix column space defined as
Γf (A,C) = Γf =
[
C⊤ · · · (CAf−1)⊤
]⊤
where f is a user-defined integer2 such that f ≥ nx. The
starting point is the following relation [18]
Yf (t) = ΓfX(t) +HfUf (t) +Vf (t) (2)
where
X(t) =
[
x(t) · · · x(t + M − 1)
]
uf (t) =
[
u⊤(t) · · · u⊤(t + f − 1)
]⊤
Uf (t) =
[
uf (t) · · · uf (t + M − 1)
]
Hf is a Toeplitz matrix composed of the system Markov
parameters, M is defined in a way compatible with the full
number of I/O measurements N and Yf (t) and Vf (t) have
the same structure as Uf (t).
B. Identification procedure
The first step3 of the propagator method consists in
applying an orthogonal projection of the row space of Yf
onto the complement of the row space of Uf
YfΠ⊥Uf = ΓfXΠ
⊥
Uf (3)
with4 Π⊥Uf = Inu −U
⊤
f (UfU⊤f )−1Uf in order to remove
the forced response and solve the unknown matrix Hf
problem.
The next step is the observability subspace estimation. It
is based on the assumptions that the system is observable
and nx is known. Assuming that (A,C) is observable, Γf
has got, at least, nx linearly independent rows. Knowing a
nx linearly independent rows set, the propagator basic idea
is to use the matrix composed of these rows as a similarity
transformation in order to fix the form of the estimated state-
space realization. The main problem is the selection of these
nx rows. When MISO observable systems are considered, the
first nx rows of Γf are linearly independent and the choice
is obvious. In the MIMO case, the problem is slightly more
difficult. Indeed, when multi-output systems are considered,
there is no way to ensure that the first nx rows of Γf are
linearly independent or that the entire system dynamics are
observable from one particular output. To get round this
difficulty, the solutions developed until now were related to
the Kronecker indices [25]. Because the state-space matrices
2When nx is known, f can be chosen equal to nx. Hereafter, equations
are described using f to tackle a more general case.
3This equation is written assuming that v = 0. The problem of noise
treatment will be explained at the end of this paragraph.
4This projection can be computed in a stable and efficient way by
resorting to the RQ factorisation [18].
6502
II. PROBLEM FORMULATION
Consider the following LTI state-space model
x(t + 1) = Ax(t) +Bu(t) (1a)
y(t) = Cx(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) are the system matrices relatively
to a certain coordinate state-space basis. It is assumed that
{x(t)}, {u(t)}, {y(t)} and {v(t)}, t ∈ Z, are all ergodic
and (weakly) stationary stochastic processes [1]. The system
order is assumed1 to be known a priori.
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 , estimate
the matrices (A,B,C) represented in a particular state-space
form and characterize the estimated parameters precision.
We start by making the following assumptions: i) the
input process {u(t)} is an ergodic and (weakly) stationary
process that is persistently exciting, ii) the output noise
{v(t)} is a zero mean noise with a finite variance statistically
uncorrelated with the input {u(t)}, iii) the system is stable
and minimal. To reach the goal mentioned above, two linked
problems are solved. Firstly, the system matrices (1) are
identified using a particular subspace-based algorithm (see
Section III). Then, the uncertainty domains of the estimated
parameters are characterized through the analysis of partic-
ular iso-level curves and an error bounded approach (see
Section IV).
III. THE PROPAGATOR METHOD
A. Basic idea of the method
Subspace-based model identification (SMI) is a particular
field of experimental modeling which has reached an inter-
esting maturity from now on [19], [16], [20], [21], [17], [18].
SMI methods are attractive because state-space realizations
can be estimated directly from I/O data without using non-
linear optimization but by applying robust linear algebra
tools such as the RQ factorization and the singular value
decomposition (SVD) [22]. Most of the SMI algorithms
lead to minimal fully-parametrized state-space realizations
[16], [17], [18]. Interesting from a numerical point of view,
the use of fully-parametrized model structure clearly over-
parametrizes the model and results in identifiability loss [23,
Chapter 2]. Furthermore, by applying a SVD, the person
using the classic subspace algorithms has difficulties to fix
the basis in which the state-space models are estimated. In
fact, this choice is realized during the reduction step (when
the system order is estimated) but the reproducibility of
the basis is not ensured when such a tool is employed. In
1This assumption is not really strong because many algorithms, mainly
based on SVD and information criteria [17], are now available to get a
reliable estimate of this parameter. Thus, for LTI systems, the system order
can be estimated beforehand.
this paper, a solution is proposed to give access to models
described in a particular structured form where the number
of elements in the state-space matrices is nx(nu + ny), i.e.
the minimal number of parameters of any MIMO system
parametrization [24].
The key problem of this identification method is the
consistent estimation, from measured I/O samples, of the
extended observability matrix column space defined as
Γf (A,C) = Γf =
[
C⊤ · · · (CAf−1)⊤
]⊤
where f is a user-defined integer2 such that f ≥ nx. The
starting point is the following relation [18]
Yf (t) = ΓfX(t) +HfUf (t) +Vf (t) (2)
where
X(t) =
[
x(t) · · · x(t + M − 1)
]
uf (t) =
[
u⊤(t) · · · u⊤(t + f − 1)
]⊤
Uf (t) =
[
uf (t) · · · uf (t + M − 1)
]
Hf is a Toeplitz matrix composed of the system Markov
parameters, M is defined in a way compatible with the full
number of I/O measurements N and Yf (t) and Vf (t) have
the same structure as Uf (t).
B. Identification procedure
The first step3 of the propagator method consists in
applying an orthogonal projection of the row space of Yf
onto the complement of the row space of Uf
YfΠ⊥Uf = ΓfXΠ
⊥
Uf (3)
with4 Π⊥Uf = Inu −U
⊤
f (UfU⊤f )−1Uf in order to remove
the forced response and solve the unknown matrix Hf
problem.
The next step is the observability subspace estimation. It
is based on the assumptions that the system is observable
and nx is known. Assuming that (A,C) is observable, Γf
has got, at least, nx linearly independent rows. Knowing a
nx linearly independent rows set, the propagator basic idea
is to use the matrix composed of these rows as a similarity
transformation in order to fix the form of the estimated state-
space realization. The main problem is the selection of these
nx rows. When MISO observable systems are considered, the
first nx rows of Γf are linearly independent and the choice
is obvious. In the MIMO case, the problem is slightly more
difficult. Indeed, when multi-output systems are considered,
there is no way to ensure that the first nx rows of Γf are
linearly independent or that the entire system dynamics are
observable from one particular output. To get round this
difficulty, the solutions developed until now were related to
the Kronecker indices [25]. Because the state-space matrices
2When nx is known, f can be chosen equal to nx. Hereafter, equations
are described using f to tackle a more general case.
3This equation is written assuming that v = 0. The problem of noise
treatment will be explained at the end of this paragraph.
4This projection can be computed in a stable and efficient way by
resorting to the RQ factorisation [18].
6502
Page 3
are unknown in our case, the determination of these indices is
problematic. In [15], it is proposed to reorganize the system
state-space representation such that a full rank matrix can
be easily extracted from the extended observability matrix.
More particularly, it is based on the following observation
(see [26, Lemma 1] for a proof)
Proposition III.1 Assume that the system is observable
and the state matrix A is non-derogatory [22]. Let κ =
[
κ1 κ2 · · · κny
]
∈ R1×ny be a vector generated ran-
domly from a uniform distribution. Then it holds with prob-
ability one that rank
{
Γnx
(
A,κ⊤C
)}
= nx.
This proposition states that the system dynamics are observ-
able from an auxiliary output ya(t) =
∑ny
i=1 κiyi(t) where yi
is the ith system output. Practically, when κj , j ∈ [1, ny], are
randomly selected, a particular matrix K defined as follows
K =
κ1 κ2 · · · κny
0 1 · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
0 0 · · · 1
∈ Rny×ny
can be introduced to substitute, e.g., ya for the first output
of the system and to obtain y¯(t) = Ky(t). Then, using this
transformation, Eq. (3) becomes
Y¯fΠ⊥Uf = Γ¯fXΠ
⊥
Uf (4)
with
Γ¯f = Γf
(
A, C¯
)
C¯ =
c¯1
c¯2
.
.
.
c¯ny
=
∑ny
j=1 κjcj
c2
.
.
.
cny
where cj , j ∈ [1, ny], are the rows of C. Noting that Γ¯f
has full column rank, a permutation matrix S (see [15] for
its construction) can be introduced in order to reorder the
rows of Γ¯f and ensure that its first nx rows are linearly
independent
SΓ¯f =
Γnx (A, c¯1)
Γf−nx (A, c¯1)Anx+1
Γnx (A, c2)
.
.
.
Γnx
(
A, cny
)
=
[
Γnx (A, c¯1)
Γcnx (A, c¯1)
]
.
By construction, rank{Γnx (A, c¯1)} = nx. Hence,
Γnx (A, c¯1) can be used as a similarity transformation T.
Furthermore, it is obvious that the complementary part
Γcnx (A, c¯1) can be written as a linear combination of
Γnx (A, c¯1). Thus a unique operator P ∈ Rnyf−nx×nx
exists, named the propagator [14], such that
Γcnx = PΓnx .
Eq. (4) becomes
SY¯fΠ⊥Uf = SΓ¯fXΠ
⊥
Uf
=
[
Inx
P
]
Γnx (A, c¯1)XΠ⊥Uf
=
[
Inx
P
]
X˜Π⊥Uf
with X˜ = TX. This relation shows that the observability
subspace can be described in a particular basis if the prop-
agator can be estimated from I/O signals. Indeed, assuming
that P is known, it holds
Υf = S⊤
[
Inx
P
]
=
C
CA
.
.
.
CAf−1
with C = CT−1 and A = TAT−1. Concerning the
estimation of P, let us introduce Z as follows
Z = SY¯fΠ⊥Uf =
[
Z1
Z2
]
=
[
Inx
P
]
X˜Π⊥Uf .
Then, under the assumptions considered until now, it is
straightforward to see that a consistent estimate of P can be
obtained by minimizing5 ‖Z2 −PZ1‖2F . Knowing Pˆ, Υˆf
can be built easily. Then, by straightforward calculations and
using the Cayley Hamilton theorem [22], it is possible to
show that
A =
[
Inx
P
]
(2 : nx + 1, :) (5)
=
0 1 0 · · · 0
0 0 1 · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 · · · 1
−a0 −a1 −a2 · · · −anx−1
(6)
c1 =
[
1 0 · · · 0
] (7)
cj =
[
Inx
P
]
((j − 1)f + 1, :) for j ∈ [2, ny] . (8)
where aj , j ∈ [0, nx − 1], are the coefficients of the char-
acteristic polynomial of A. Once A and C are known, the
matrix B can be obtained by linear regression [18]. Observ-
ing A, B and C, the number of parameters6 is nx(nu +ny),
i.e. the minimal number of parameters of any MIMO system
parametrization [24]. So, thanks to the minimal state-space
matrices structure provided by the propagator method, the
uncertainty areas description will be easier.
In the noisy data case, this method can be adapted by
introducing an instrumental variable, as proposed, e.g., in
the MOESP approach [18]. The choice of the instruments
depends on the inputs and disturbances properties (see [18]
for details).
5‖.‖F is the Frobenius norm.
6Notice that B is a fully-parametrized matrix.
6503
problematic. In [15], it is proposed to reorganize the system
state-space representation such that a full rank matrix can
be easily extracted from the extended observability matrix.
More particularly, it is based on the following observation
(see [26, Lemma 1] for a proof)
Proposition III.1 Assume that the system is observable
and the state matrix A is non-derogatory [22]. Let κ =
[
κ1 κ2 · · · κny
]
∈ R1×ny be a vector generated ran-
domly from a uniform distribution. Then it holds with prob-
ability one that rank
{
Γnx
(
A,κ⊤C
)}
= nx.
This proposition states that the system dynamics are observ-
able from an auxiliary output ya(t) =
∑ny
i=1 κiyi(t) where yi
is the ith system output. Practically, when κj , j ∈ [1, ny], are
randomly selected, a particular matrix K defined as follows
K =
κ1 κ2 · · · κny
0 1 · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
0 0 · · · 1
∈ Rny×ny
can be introduced to substitute, e.g., ya for the first output
of the system and to obtain y¯(t) = Ky(t). Then, using this
transformation, Eq. (3) becomes
Y¯fΠ⊥Uf = Γ¯fXΠ
⊥
Uf (4)
with
Γ¯f = Γf
(
A, C¯
)
C¯ =
c¯1
c¯2
.
.
.
c¯ny
=
∑ny
j=1 κjcj
c2
.
.
.
cny
where cj , j ∈ [1, ny], are the rows of C. Noting that Γ¯f
has full column rank, a permutation matrix S (see [15] for
its construction) can be introduced in order to reorder the
rows of Γ¯f and ensure that its first nx rows are linearly
independent
SΓ¯f =
Γnx (A, c¯1)
Γf−nx (A, c¯1)Anx+1
Γnx (A, c2)
.
.
.
Γnx
(
A, cny
)
=
[
Γnx (A, c¯1)
Γcnx (A, c¯1)
]
.
By construction, rank{Γnx (A, c¯1)} = nx. Hence,
Γnx (A, c¯1) can be used as a similarity transformation T.
Furthermore, it is obvious that the complementary part
Γcnx (A, c¯1) can be written as a linear combination of
Γnx (A, c¯1). Thus a unique operator P ∈ Rnyf−nx×nx
exists, named the propagator [14], such that
Γcnx = PΓnx .
Eq. (4) becomes
SY¯fΠ⊥Uf = SΓ¯fXΠ
⊥
Uf
=
[
Inx
P
]
Γnx (A, c¯1)XΠ⊥Uf
=
[
Inx
P
]
X˜Π⊥Uf
with X˜ = TX. This relation shows that the observability
subspace can be described in a particular basis if the prop-
agator can be estimated from I/O signals. Indeed, assuming
that P is known, it holds
Υf = S⊤
[
Inx
P
]
=
C
CA
.
.
.
CAf−1
with C = CT−1 and A = TAT−1. Concerning the
estimation of P, let us introduce Z as follows
Z = SY¯fΠ⊥Uf =
[
Z1
Z2
]
=
[
Inx
P
]
X˜Π⊥Uf .
Then, under the assumptions considered until now, it is
straightforward to see that a consistent estimate of P can be
obtained by minimizing5 ‖Z2 −PZ1‖2F . Knowing Pˆ, Υˆf
can be built easily. Then, by straightforward calculations and
using the Cayley Hamilton theorem [22], it is possible to
show that
A =
[
Inx
P
]
(2 : nx + 1, :) (5)
=
0 1 0 · · · 0
0 0 1 · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 · · · 1
−a0 −a1 −a2 · · · −anx−1
(6)
c1 =
[
1 0 · · · 0
] (7)
cj =
[
Inx
P
]
((j − 1)f + 1, :) for j ∈ [2, ny] . (8)
where aj , j ∈ [0, nx − 1], are the coefficients of the char-
acteristic polynomial of A. Once A and C are known, the
matrix B can be obtained by linear regression [18]. Observ-
ing A, B and C, the number of parameters6 is nx(nu +ny),
i.e. the minimal number of parameters of any MIMO system
parametrization [24]. So, thanks to the minimal state-space
matrices structure provided by the propagator method, the
uncertainty areas description will be easier.
In the noisy data case, this method can be adapted by
introducing an instrumental variable, as proposed, e.g., in
the MOESP approach [18]. The choice of the instruments
depends on the inputs and disturbances properties (see [18]
for details).
5‖.‖F is the Frobenius norm.
6Notice that B is a fully-parametrized matrix.
6503
Page 4
IV. UNCERTAINTY DOMAIN DETERMINATION
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
description of the coefficients ai, bji and cki , i ∈ [0, nx − 1],
j ∈ [1, nu], k ∈ [2, ny] (the state-space matrices A, B
and C parameters) is performed. The developed approach
basic idea is to study particular ellipsoidal iso-level curves
by adapting a bounded-error method proposed in [13]. This
original method has been developed assuming that the system
model is described via a linear regression. Because this
assumption is not satisfied with the estimated state-space
model, it is firstly necessary to transform the estimated state-
space form into a representation linear in the parameters.
Notice right now that the parameters of this modified model
are combinations of the initial ai, bji and cki , i ∈ [0, nx − 1],
j ∈ [1, nu], k ∈ [2, ny]. Then, the method developed in [13]
is adapted to this MIMO representation linear in the param-
eters. Finally, the uncertainty domain of the coefficients ai,
bji and cki are deduced from the descriptions obtained by the
I/O model. All these phases are presented in the following.
A. Input-output model description
1) General case: In this paragraph, a new representation
of the estimated state-space model is proposed. More partic-
ularly, an input-output (I/O) ARMAX-like model is obtained
by getting rid of the state in Eq. (1). All the calculations are
available in [27]. Due to the lack of space, only the main
proposition is given hereafter. More precisely,
Proposition IV.1 Consider the state-space representation
(1). Then, it can be described by an ARMAX-like model
y(t) = −
(
a⊤ ⊗ Iny
)
ynx(t − nx)
+ Funx(t − nx) + η(t) = Φ⊤(t)θ + η(t) (9)
with
Φ⊤(t) = [−y(t−nx) ··· −y(t−1) unx(t−nx)⊗Iny ] (10)
θ =
[
a
vec(F)
]
(11)
a⊤ =
[
a0 a1 · · · anx−1
] (12)
F =
(
a⊤ ⊗ Iny
)
Hnx +C∆nx ∈ Rny×nxnu (13)
∆nx =
[
Anx−1B · · · AB B
]
∈ Rnx×nxnu (14)
and η(t) a noise corresponding to a filtered version of v.
Furthermore, the order of the I/O model (9) is equal to nx if
and only if the initial system (1) is minimal and the matrix
A is non-derogatory.
It is important to note that, in the following, condition “A is
non-derogatory” will be satisfied. Indeed, this transformation
will be applied to the estimated model, i.e. with a matrix A
in a companion form (see Eq. (6)).
2) Case study: Consider the following SISO state-space
model
x(t + 1) =
0 1 0
0 0 1
−a0 −a1 −a2
x(t) +
b10
b11
b12
u(t) (15a)
y(t) =
[
1 0 0
]
x(t) (15b)
Then, the I/O model satisfies
y(t) = −a2y(t−1)−a1y(t−2)−a0y(t−3)+ b10u(t−1)
+
(
a2b10 + b11
)
u(t− 2) +
(
a1b10 + a2b11 + b12
)
u(t− 3).
It is obvious that this I/O model is linear in θ but non-
linear in ai, bji and cki . Remember that the final goal is to
get the uncertainty domains description of the state-space
representation coefficients.
B. Uncertainty domain description for a model linear in the
parameters
In this paragraph, the uncertainty domains description
of the parameters θi composing θ is considered7. Looking
closes at Eq. (9), it is obvious that the least-squares estima-
tion of this parameter vector can be obtained by minimizing
the following cost function [1, Appendix II]
J(θ¯) = 12
(
yM −ΨM θ¯
)⊤ (
yM −ΨM θ¯
)
where θ¯ is an arbitrary estimate of θ and with
yM =
y(1)
.
.
.
y(M)
ΨM =
Φ⊤(1)
.
.
.
Φ⊤(M)
.
Then, assuming that RΨ = Ψ⊤MΨM is invertible, it can be
shown that
J(θ¯) = 12y
⊤
M
(
I−ΨM
(
Ψ⊤MΨM
)−1Ψ⊤M
)
yM
+ 12
(
θ¯ − θls
)⊤Ψ⊤MΨM
(
θ¯ − θls
)
= Jmin +
1
2dθ¯
⊤RΨdθ¯ (16)
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 rhs of Eq. (16) is the minimum value.
Furthermore, if we introduce V (θ¯) = J(θ¯) − Jmin, it is
obvious that
V (θ¯) = 12dθ¯
⊤RΨdθ¯ (17)
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
7The link with the coefficients ai, bji and cki is postponed for §IV-C.
6504
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
description of the coefficients ai, bji and cki , i ∈ [0, nx − 1],
j ∈ [1, nu], k ∈ [2, ny] (the state-space matrices A, B
and C parameters) is performed. The developed approach
basic idea is to study particular ellipsoidal iso-level curves
by adapting a bounded-error method proposed in [13]. This
original method has been developed assuming that the system
model is described via a linear regression. Because this
assumption is not satisfied with the estimated state-space
model, it is firstly necessary to transform the estimated state-
space form into a representation linear in the parameters.
Notice right now that the parameters of this modified model
are combinations of the initial ai, bji and cki , i ∈ [0, nx − 1],
j ∈ [1, nu], k ∈ [2, ny]. Then, the method developed in [13]
is adapted to this MIMO representation linear in the param-
eters. Finally, the uncertainty domain of the coefficients ai,
bji and cki are deduced from the descriptions obtained by the
I/O model. All these phases are presented in the following.
A. Input-output model description
1) General case: In this paragraph, a new representation
of the estimated state-space model is proposed. More partic-
ularly, an input-output (I/O) ARMAX-like model is obtained
by getting rid of the state in Eq. (1). All the calculations are
available in [27]. Due to the lack of space, only the main
proposition is given hereafter. More precisely,
Proposition IV.1 Consider the state-space representation
(1). Then, it can be described by an ARMAX-like model
y(t) = −
(
a⊤ ⊗ Iny
)
ynx(t − nx)
+ Funx(t − nx) + η(t) = Φ⊤(t)θ + η(t) (9)
with
Φ⊤(t) = [−y(t−nx) ··· −y(t−1) unx(t−nx)⊗Iny ] (10)
θ =
[
a
vec(F)
]
(11)
a⊤ =
[
a0 a1 · · · anx−1
] (12)
F =
(
a⊤ ⊗ Iny
)
Hnx +C∆nx ∈ Rny×nxnu (13)
∆nx =
[
Anx−1B · · · AB B
]
∈ Rnx×nxnu (14)
and η(t) a noise corresponding to a filtered version of v.
Furthermore, the order of the I/O model (9) is equal to nx if
and only if the initial system (1) is minimal and the matrix
A is non-derogatory.
It is important to note that, in the following, condition “A is
non-derogatory” will be satisfied. Indeed, this transformation
will be applied to the estimated model, i.e. with a matrix A
in a companion form (see Eq. (6)).
2) Case study: Consider the following SISO state-space
model
x(t + 1) =
0 1 0
0 0 1
−a0 −a1 −a2
x(t) +
b10
b11
b12
u(t) (15a)
y(t) =
[
1 0 0
]
x(t) (15b)
Then, the I/O model satisfies
y(t) = −a2y(t−1)−a1y(t−2)−a0y(t−3)+ b10u(t−1)
+
(
a2b10 + b11
)
u(t− 2) +
(
a1b10 + a2b11 + b12
)
u(t− 3).
It is obvious that this I/O model is linear in θ but non-
linear in ai, bji and cki . Remember that the final goal is to
get the uncertainty domains description of the state-space
representation coefficients.
B. Uncertainty domain description for a model linear in the
parameters
In this paragraph, the uncertainty domains description
of the parameters θi composing θ is considered7. Looking
closes at Eq. (9), it is obvious that the least-squares estima-
tion of this parameter vector can be obtained by minimizing
the following cost function [1, Appendix II]
J(θ¯) = 12
(
yM −ΨM θ¯
)⊤ (
yM −ΨM θ¯
)
where θ¯ is an arbitrary estimate of θ and with
yM =
y(1)
.
.
.
y(M)
ΨM =
Φ⊤(1)
.
.
.
Φ⊤(M)
.
Then, assuming that RΨ = Ψ⊤MΨM is invertible, it can be
shown that
J(θ¯) = 12y
⊤
M
(
I−ΨM
(
Ψ⊤MΨM
)−1Ψ⊤M
)
yM
+ 12
(
θ¯ − θls
)⊤Ψ⊤MΨM
(
θ¯ − θls
)
= Jmin +
1
2dθ¯
⊤RΨdθ¯ (16)
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 rhs of Eq. (16) is the minimum value.
Furthermore, if we introduce V (θ¯) = J(θ¯) − Jmin, it is
obvious that
V (θ¯) = 12dθ¯
⊤RΨdθ¯ (17)
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
7The link with the coefficients ai, bji and cki is postponed for §IV-C.
6504
Page 5
too conservative solution. This level is obtained hereafter
by using an analogy with the stochastic framework. More
precisely, under the Gaussian case 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 in Eq. (17).
This analogy leads to choose the level JD as follows
JD − Jmin = n2σ2.
Unfortunately, this relation depends on the prior knowledge
of σ2. In order to circumvent this difficulty, a bounded-
error approach is used8. 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 ε = y − yˆ(θ¯) are
the residuals. Thus, a good noise effect approximation can
be deduced from ε and we fix ℓ as follows
ℓ = max {|ε(t)|} .
C. Extension to a model non-linear in the parameters
The final goal of the developed method is to describe the
uncertainty domains of the coefficients ai, bji and cki , i ∈
[0, nx − 1], j ∈ [1, nu], k ∈ [2, ny]. To reach this goal, let
us introduce θss as the parameters vector of the estimated
state-space form, i.e.
θss =
[
−a0 · · · −anx−1 b10 · · · bnu0 b11 · · ·
bnunx−1 c
2
0 · · · c
ny
0 c21 · · · c
ny
nx−1
]
.
Then, we know that θss and θ can be related via a mapping
f known a priori (see Eq. (11)), i.e. θ = f(θss). Using a
Taylor series expansion, the following first order approxima-
tion is satisfied
dθ ≈
(
∂f
∂θss
)
θss=θˆss
dθss
where θˆss is the estimated state-space parameters vector.
Combining this relation with Eq. (16), it holds
J ≈ Jmin + dθ⊤ssRdθss
with
R =
(
∂f
∂θss
)⊤
θss=θˆss
RΨ
(
∂f
∂θss
)
θss=θˆss
.
8In this method, the only assumption is to have bounded residuals.
The uncertainty domains are again paraboloids whose main
directions are given by R, i.e. a modified version of RΨ.
Hence, the procedure explained in the LP case can be easily
extended for θss.
D. Case study
In order to illustrate the developed approach, assume that
the identified model is described by the state-space form (15).
Then,
θss =
[
−a0 −a1 −a2 b10 b11 b12
]
θ =
[
a0 a1 a2 a1b10 + a2b11 + b12 a2b10 + b11 b10
]
.
Furthermore,
(
∂f
∂θss
)
=
−1 0 0 0 0 0
0 −1 0 0 0 0
0 0 −1 0 0 0
0 −b10 −b11 a1 a2 1
0 0 −b10 a2 1 0
0 0 0 1 0 0
.
V. SIMULATION EXAMPLE
In order to show the performances 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 =
2
4
0 1 0
0 0 1
0.1 0.2 −0.1
3
5B =
2
4
7.3 0.2
−0.3 −6.9
−14.6 −2.0
3
5C =
»
1.0 0 0
1.0 −0.3 0.6
–
.
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 equals
20dB.
The propagator method is applied to estimate the system
matrices. To reduce the noise effect, past inputs are used as
instruments (as for the PI MOESP algorithm [18]). The past
and future horizons are chosen equal to 5.
Identification results are plotted in Figure 1. Firstly, these
two graphs show that the estimates obtained using the
propagator method are accurate. Indeed, in each plot, the
black cross (+), symbolizing the mean value of the esti-
mated parameters, almost matches the red one (×) which
corresponds to the real parameters9.
As far as the uncertainty domains are concerned, the level
ℓ is estimated such that ℓ = 1.11. 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. To assess
the quality of these uncertainty domains and to get rid of
the drawing of all the uncertainty domains, a failure rate
measure, defined as the percentage of realizations for which
the system parameters are outside of the ellipsoid D centered
in the estimated parameters vector, is used (see black disc
9Knowing the state-space matrices, it is straightforward to compute the
similarity transformation Γnx (A, c¯1) and to transform the system state-
space form into the structure (6)-(8).
6505
by using an analogy with the stochastic framework. More
precisely, under the Gaussian case 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 in Eq. (17).
This analogy leads to choose the level JD as follows
JD − Jmin = n2σ2.
Unfortunately, this relation depends on the prior knowledge
of σ2. In order to circumvent this difficulty, a bounded-
error approach is used8. 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 ε = y − yˆ(θ¯) are
the residuals. Thus, a good noise effect approximation can
be deduced from ε and we fix ℓ as follows
ℓ = max {|ε(t)|} .
C. Extension to a model non-linear in the parameters
The final goal of the developed method is to describe the
uncertainty domains of the coefficients ai, bji and cki , i ∈
[0, nx − 1], j ∈ [1, nu], k ∈ [2, ny]. To reach this goal, let
us introduce θss as the parameters vector of the estimated
state-space form, i.e.
θss =
[
−a0 · · · −anx−1 b10 · · · bnu0 b11 · · ·
bnunx−1 c
2
0 · · · c
ny
0 c21 · · · c
ny
nx−1
]
.
Then, we know that θss and θ can be related via a mapping
f known a priori (see Eq. (11)), i.e. θ = f(θss). Using a
Taylor series expansion, the following first order approxima-
tion is satisfied
dθ ≈
(
∂f
∂θss
)
θss=θˆss
dθss
where θˆss is the estimated state-space parameters vector.
Combining this relation with Eq. (16), it holds
J ≈ Jmin + dθ⊤ssRdθss
with
R =
(
∂f
∂θss
)⊤
θss=θˆss
RΨ
(
∂f
∂θss
)
θss=θˆss
.
8In this method, the only assumption is to have bounded residuals.
The uncertainty domains are again paraboloids whose main
directions are given by R, i.e. a modified version of RΨ.
Hence, the procedure explained in the LP case can be easily
extended for θss.
D. Case study
In order to illustrate the developed approach, assume that
the identified model is described by the state-space form (15).
Then,
θss =
[
−a0 −a1 −a2 b10 b11 b12
]
θ =
[
a0 a1 a2 a1b10 + a2b11 + b12 a2b10 + b11 b10
]
.
Furthermore,
(
∂f
∂θss
)
=
−1 0 0 0 0 0
0 −1 0 0 0 0
0 0 −1 0 0 0
0 −b10 −b11 a1 a2 1
0 0 −b10 a2 1 0
0 0 0 1 0 0
.
V. SIMULATION EXAMPLE
In order to show the performances 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 =
2
4
0 1 0
0 0 1
0.1 0.2 −0.1
3
5B =
2
4
7.3 0.2
−0.3 −6.9
−14.6 −2.0
3
5C =
»
1.0 0 0
1.0 −0.3 0.6
–
.
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 equals
20dB.
The propagator method is applied to estimate the system
matrices. To reduce the noise effect, past inputs are used as
instruments (as for the PI MOESP algorithm [18]). The past
and future horizons are chosen equal to 5.
Identification results are plotted in Figure 1. Firstly, these
two graphs show that the estimates obtained using the
propagator method are accurate. Indeed, in each plot, the
black cross (+), symbolizing the mean value of the esti-
mated parameters, almost matches the red one (×) which
corresponds to the real parameters9.
As far as the uncertainty domains are concerned, the level
ℓ is estimated such that ℓ = 1.11. 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. To assess
the quality of these uncertainty domains and to get rid of
the drawing of all the uncertainty domains, a failure rate
measure, defined as the percentage of realizations for which
the system parameters are outside of the ellipsoid D centered
in the estimated parameters vector, is used (see black disc
9Knowing the state-space matrices, it is straightforward to compute the
similarity transformation Γnx (A, c¯1) and to transform the system state-
space form into the structure (6)-(8).
6505
Page 6
(•) in Fig. 1). The failure ratio equals 1.7 % for the couple
(−a0,−a1) and 2.4 % for (−a2, b10). These values show that
the way the threshold is fixed leads to reliable uncertainty
domains.
0.161 0.162 0.163 0.164 0.165 0.166 0.167
0.237
0.238
0.239
0.24
0.241
0.242
0.243
−a0
−
a 1
−0.103 −0.102 −0.101 −0.1 −0.099 −0.098 −0.0977.33
7.34
7.35
7.36
7.37
7.38
7.39
7.4
7.41
7.42
7.43
−a2
b1 0
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.
VI. CONCLUSIONS
In this paper, the problem of uncertainty domain deter-
mination for MIMO LTI state-space systems is considered.
To solve this problem, a particular subspace-based method is
used. This technique leads to an estimated state-space form
with fixed structure and a minimal number of parameters.
The description of the uncertainty domains is then realized
in three steps. Firstly, the estimated model is transformed into
an ARMAX-like representation. The uncertainty domains
of the I/O model parameters are then described using a
bounded-error approach. The uncertainty domains of the
initial model parameters are finally deduced from those
obtained in the former step by adapting particular cost
functions. The experimental results have emphasized the
reliability of the developed method.
REFERENCES
[1] L. Ljung, System identification. Theory for the user, 2nd ed. Upper
Saddle River: Prentice Hall, 1999.
[2] 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, 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] X. Bombois, “Connecting prediction error identification and robust
control analysis: a new framework,” Ph.D. dissertation, Université
Catholique de Louvain, Louvain la Neuve, Belgium, 2000.
[6] R. Hakvoort, “System identification for robust process control -
nominal models and error bounds,” Ph.D. dissertation, Delft University
of Technology, Delft, The Netherlands, 1994.
[7] D. de Vries, “Identification of model uncertainty for control design,”
Ph.D. dissertation, Delft University of Technology, Delft, The Nether-
lands, 1994.
[8] L. Giarré and M. Milanese, “Model quality evaluation in H2 identifica-
tion,” IEEE Transactions on Automatic Control, vol. 42, pp. 691–698,
1997.
[9] L. Giarré, M. Milanese, and M. Taragna, “H∞ identification and
model quality evaluation,” IEEE Transactions on Automatic Control,
vol. 42, pp. 188–199, 1997.
[10] E. Walter, J. Norton, H. Piet-Lahanier, and M. Milanese, Bounding
Approaches to System Identification. Perseus Publishing, 1996.
[11] 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.
[12] 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.
[13] C. Baron, T. Poinot, and J. 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, 2001.
[14] J. Munier and G. Delisle, “Spatial analysis using new properties of
the cross spectral matrix,” IEEE Transactions on Signal Processing,
vol. 39, pp. 746–749, 1991.
[15] 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.
[16] P. Van Overschee and B. De Moor, Subspace identification for linear
systems. theory, implementation, applications. Kluwer Academic
Publishers, 1996.
[17] T. Katayama, Subspace methods for system identification. Springer
Verlag, 2005.
[18] M. Verhaegen and V. Verdult, Filtering and system identification: a
least squares approach. Cambridge University Press, 2007.
[19] S. Kukreja, B. Haverkamp, D. Westwick, R. Kearney, H. Galiana, and
M. Verhaegen., “Subspace identification method for ankle mechanics,”
IEEE Engineering in Medecine and Biology Society, vol. 17, pp. 1413–
1414, 1995.
[20] M. Abdelghani, M. Verhaegen, P. Van Overschee, and B. De Moor,
“Comparison study of subspace identification methods applied to flex-
ible structures,” Mechanical Systems and Signal Processing, vol. 12,
pp. 679–692, 1998.
[21] W. Favoreel, B. De Moor, and P. Van Overschee, “Subspace state space
system identification for industrial processes,” Journal of Process
Control, vol. 10, pp. 149–155, 2000.
[22] R. Horn and C. Johnson, Matrix analysis. Cambridge University
Press, 1990.
[23] T. McKelvey, “Identification of state space models from time and
frequency data,” Ph.D. dissertation, Linköping University, Linköping,
Sweden, 1995.
[24] T. McKelvey and A. Helmersson, “System identification using an over-
parametrized model class improving the optimization algorithm,” in
Proceedings of the 36th IEEE Conference on Decision and Control,
San Diego, California USA, December 1997.
[25] T. Kailath, Linear Systems. Engelwood Cliffs: Prentice Hall, 1980.
[26] 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.
[27] ——, “Identification of multivariable canonical state-space represen-
tations,” in Proceedings of the 15th IFAC Symposium on System
Identification, Saint Malo, France, July 2009.
6506
(−a0,−a1) and 2.4 % for (−a2, b10). These values show that
the way the threshold is fixed leads to reliable uncertainty
domains.
0.161 0.162 0.163 0.164 0.165 0.166 0.167
0.237
0.238
0.239
0.24
0.241
0.242
0.243
−a0
−
a 1
−0.103 −0.102 −0.101 −0.1 −0.099 −0.098 −0.0977.33
7.34
7.35
7.36
7.37
7.38
7.39
7.4
7.41
7.42
7.43
−a2
b1 0
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.
VI. CONCLUSIONS
In this paper, the problem of uncertainty domain deter-
mination for MIMO LTI state-space systems is considered.
To solve this problem, a particular subspace-based method is
used. This technique leads to an estimated state-space form
with fixed structure and a minimal number of parameters.
The description of the uncertainty domains is then realized
in three steps. Firstly, the estimated model is transformed into
an ARMAX-like representation. The uncertainty domains
of the I/O model parameters are then described using a
bounded-error approach. The uncertainty domains of the
initial model parameters are finally deduced from those
obtained in the former step by adapting particular cost
functions. The experimental results have emphasized the
reliability of the developed method.
REFERENCES
[1] L. Ljung, System identification. Theory for the user, 2nd ed. Upper
Saddle River: Prentice Hall, 1999.
[2] 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, 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] X. Bombois, “Connecting prediction error identification and robust
control analysis: a new framework,” Ph.D. dissertation, Université
Catholique de Louvain, Louvain la Neuve, Belgium, 2000.
[6] R. Hakvoort, “System identification for robust process control -
nominal models and error bounds,” Ph.D. dissertation, Delft University
of Technology, Delft, The Netherlands, 1994.
[7] D. de Vries, “Identification of model uncertainty for control design,”
Ph.D. dissertation, Delft University of Technology, Delft, The Nether-
lands, 1994.
[8] L. Giarré and M. Milanese, “Model quality evaluation in H2 identifica-
tion,” IEEE Transactions on Automatic Control, vol. 42, pp. 691–698,
1997.
[9] L. Giarré, M. Milanese, and M. Taragna, “H∞ identification and
model quality evaluation,” IEEE Transactions on Automatic Control,
vol. 42, pp. 188–199, 1997.
[10] E. Walter, J. Norton, H. Piet-Lahanier, and M. Milanese, Bounding
Approaches to System Identification. Perseus Publishing, 1996.
[11] 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.
[12] 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.
[13] C. Baron, T. Poinot, and J. 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, 2001.
[14] J. Munier and G. Delisle, “Spatial analysis using new properties of
the cross spectral matrix,” IEEE Transactions on Signal Processing,
vol. 39, pp. 746–749, 1991.
[15] 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.
[16] P. Van Overschee and B. De Moor, Subspace identification for linear
systems. theory, implementation, applications. Kluwer Academic
Publishers, 1996.
[17] T. Katayama, Subspace methods for system identification. Springer
Verlag, 2005.
[18] M. Verhaegen and V. Verdult, Filtering and system identification: a
least squares approach. Cambridge University Press, 2007.
[19] S. Kukreja, B. Haverkamp, D. Westwick, R. Kearney, H. Galiana, and
M. Verhaegen., “Subspace identification method for ankle mechanics,”
IEEE Engineering in Medecine and Biology Society, vol. 17, pp. 1413–
1414, 1995.
[20] M. Abdelghani, M. Verhaegen, P. Van Overschee, and B. De Moor,
“Comparison study of subspace identification methods applied to flex-
ible structures,” Mechanical Systems and Signal Processing, vol. 12,
pp. 679–692, 1998.
[21] W. Favoreel, B. De Moor, and P. Van Overschee, “Subspace state space
system identification for industrial processes,” Journal of Process
Control, vol. 10, pp. 149–155, 2000.
[22] R. Horn and C. Johnson, Matrix analysis. Cambridge University
Press, 1990.
[23] T. McKelvey, “Identification of state space models from time and
frequency data,” Ph.D. dissertation, Linköping University, Linköping,
Sweden, 1995.
[24] T. McKelvey and A. Helmersson, “System identification using an over-
parametrized model class improving the optimization algorithm,” in
Proceedings of the 36th IEEE Conference on Decision and Control,
San Diego, California USA, December 1997.
[25] T. Kailath, Linear Systems. Engelwood Cliffs: Prentice Hall, 1980.
[26] 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.
[27] ——, “Identification of multivariable canonical state-space represen-
tations,” in Proceedings of the 15th IFAC Symposium on System
Identification, Saint Malo, France, July 2009.
6506
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
2 Readers on Mendeley
by Discipline
by Academic Status
50% Student (Bachelor)
50% Associate Professor
by Country
50% Brazil
50% France


