RECURSIVE MOESP TYPE SUBSPACE IDENTIFICATION

ALGORITHM

Catarina J. M. Delgado

Faculdade de Economia do Porto

Rua Dr. Roberto Frias, 4200-464 Porto - Portugal

P. Lopes dos Santos

Faculdade de Engenharia da Universidade do Porto

Rua Dr. Roberto Frias, 4200-465 Porto - Portugal

Keywords:

System Identiﬁcation, Subspace methods, State-space models, Linear Systems, Parameter estimation.

Abstract:

In this paper two recursive algorithms, based on MOESP type subspace identiﬁcation, are presented in two

versions. The main idea was to show that we can represent subspace identiﬁcation methods as sequences of

least squares problems and implement them through sequences of modiﬁed Householder algorithms. There-

fore, it is possible to develop iterative algorithms with most of the advantages of this kind of methods, and still

improve the numerical efﬁciency, in order to deal with real-time applications and minimize the computational

burden.

1 INTRODUCTION

The Subspace Identiﬁcation algorithms aim to esti-

mate, from measured input / output data sequences

({u

} and {y

}, respectively), the system described

by:

k+1

= Ax

+ Bu

+ Ke

= Cx

+ Du

+ e

(1)

= R

> 0 (2)

where A ∈ R

n×n

, B ∈ R

n×m

, C ∈ R

l×n

, D ∈

l×m

, K ∈ R

n×l

and x

∈ R

. The sequence

} ∈ R

is a white noise stochastic process and the

input data sequence is assumed to be a persistently ex-

citing quasi-stationary deterministic sequence(Ljung,

1987).

Most of the past contributions ((Larimore, 1990),

(Moonen et al., 1989), (van Overschee and de Moor,

1996), (Verhaegen, 1994), (Viberg et al., 1991))

are centered in the off-line identiﬁcation of time-

invariant systems ﬁeld, by developing non-recursive

algorithms. Even though the on-line identiﬁcation can

be applied to a wider class of systems (time-variant,

non-linear...), it has been a little forgotten due to its

high computation and storage costs. When a particu-

lar subspace off-line identiﬁcation algorithm (such as

the algorithm introduced by Van Overschee and De

Moor (van Overschee and de Moor, 1996)) is adapted

to a recursive identiﬁcation algorithm ((Cho et al.,

1994), (Ljung, 1983), (Gustafsson et al., 1998), (Oku

and Kimura, 2000)), not only the computational bur-

den is reduced, but also the application scope is en-

larged, preserving most of the advantages of this kind

of methods.

This paper is organized as follows: in section 2,

the notation is introduced and, in section 3, two ba-

sic MOESP-type algorithms are described: the PO-

MOESP (Verhaegen, 1994) and the R-MOESP (van

Overschee and de Moor, 1996). Sections 4 and 5 are

dedicated to the main results: section 4 introduces

the 2QR recursive version for the MOESP algorithms

and section 5 the 1QR version. Finally, in section 6,

some simulation results are shown and, in section 7,

the conclusions are presented.

2 NOTATION

We deﬁne the following block Hankel matrices, built

with input data:





(0)

...

(i−1)





, U





(0)

...

(i)









(i)

...

(2i−1)





, U

f−





(i+1)

...

(2i)





where U

(i)

= [

... u

i+j−1

] ∈ R

m×j

and the

subscripts p and f denote past and future data, respec-

176

Delgado C. and dos Santos P. (2004).

RECURSIVE MOESP TYPE SUBSPACE IDENTIFICATION ALGORITHM.

In Proceedings of the First International Conference on Informatics in Control, Automation and Robotics, pages 176-182

DOI: 10.5220/0001146801760182

 SciTePress

tively. In a similar way, Y

(i)

∈ R

l×j

, Y

and

f−

are block Hankel matrices, built with output data.

With both input and output data we can also deﬁne

the block hankel matrices

, H

HU =

, H

−

f−

and the orthogonal projection matrices

= Y

i+1

= Y

−

= Y

−

where Π

= M

M = M

M denotes

the orthogonal projection into the row space of M and

(.)

denotes the Moore-Penrose pseudoinverse

Other important matrices are the extended observ-

ability matrix

∈ R

li×n

and the block Toeplitz

matrix

∈ R

li×mi

, built with Markov parame-

ters (Delgado and dos Santos, 2004b).

3 MOESP-TYPE ALGORITHMS

The MOESP algorithms can be described as se-

quences of two main steps (van Overschee and

de Moor, 1996). First, the model order n and a ba-

sis for the column space of Γ

are determined:

1/2

, where U

and S

are given by the Singular

Value Decomposition of Z

⊥

= Z

− Z

⊥

= USV

= (3)

= [

⊥

]

0 0

0 S

⊥





⊥





= U

+ U

⊥

≈

≈ U

When the measurements are noise corrupted, it may

not be straightforward to distinguish the ”nonzero”

singular values (S

) from the ”zero” singular values

⊥

). Therefore, a better estimate of the system or-

der is obtained by trying different values and then

comparing the simulation errors (van Overschee and

de Moor, 1996).

In the second main step, the system matrices are es-

timated. There are two approaches to estimate those

matrices: the method introduced by Van Overschee,

with the R-MOESP algorithm (van Overschee and

de Moor, 1996), and the method introduced by Ver-

haegen, with the PO-MOESP algorithm (Verhaegen,

1994).

3.1 R-MOESP Algorithm.

When we know an estimate of Γ

, matrices A and C

are obtained by solving the following linear equation,

in a least squares sense (van Overschee and de Moor,

1996):

i−1

i+1

(i)

¸·

+ ρ (4)

where Γ

i−1

l(i−1)

. The matrix K

then used to estimate B and D, since

= (5)

B Γ

i−1

D 0

−

... N

¸ ¸

Lopes dos Santos and Martins de Carvalho (dos

Santos and de Carvalho, 2003) have shown that

can be written as K

= K

−A

−1

. Therefore, we can work

only with:

C(B,D)

= (6)

... N

¸ ¸

where

= [

− L

−L

... −L

] G

= [

−L

... −L

] G

(k > 1)

= [

... LA

] = AΓ

= [

... L

] = CΓ

0 Γ

i−1

Equation (6) can be rewritten as:

...

(7)

and B and D estimated in the least squares sense.

3.2 PO-MOESP Algorithm

Another approach to determine matrices A and C is

to use the shift-invariance structure in Γ

(Verhaegen,

1994):

C = [

]

(8)

A = Γ

i−1

(9)

RECURSIVE MOESP TYPE SUBSPACE IDENTIFICATION ALGORITHM

177

where Γ

0 I

l(i−1)

. These estimates of A

and C are the same as the estimates produced by the

method of Van Overschee and De Moor ((Chiuso and

Picci, 2001), (Delgado and dos Santos, 2004a), (dos

Santos and de Carvalho, 2004)).

Then, matrices B an D can be estimated by linear

regression:

vec(by

) = [M

]

(10)

where

= (u

⊗ I

) (11)

(k−1)

t=1

⊗ CA

k−t−1

) (12)

= vec(D) (13)

= vec(B) (14)

and by

is the simulated output, given by

= Du

k−1

t=1

k−t−1

= by

+ by

(15)

with

vec(by

) = vec(Du

) = (16)

= (u

⊗ I

)vec(D)

vec(by

) = vec

k−1

t=1

k−t−1

(17)

k−1

t=1

⊗ CA

k−t−1

vec(B)

4 RECURSIVE 2QR APPROACH

4.1 The weighted projection

The estimation of Γ

from the singular value decom-

position of Z

⊥

is implemented through the LQ

decomposition of

= LQ

= (18)

0 0

where L

∈ R

mi×mi

, L

∈ R

(m+l)i×mi

∈ R

(m+l)i×(m+l)i

and Q

∈ R

mi×j

, Q

∈

(m+l)i×j

, Q

∈ R

(j−(2m+l)i)×j

. Matrices L

and L

are two lower triangular matrices and ma-

trix Q

is orthonormal and

= I

⇔ Q

I, if p = q

0, if p 6= q

(19)

k=1

= I

(20)

From this decomposition,

= U

)

−1

= Q

(21)

(22)

⊥

= I

− Π

(23)

⊥

= I

− Π

= Q

(24)

If we consider the least squares problem

= [θ

]

(25)

we can rewrite it as

= [θ

]

0 0

(26)

where B

⊥

= Y

Then,

= Y

= [B

]

(27)

= B

+ B

As to the weighted projection, it can be expressed as

⊥

= θ

⊥

= (28)

= [

]

⊥

= Z

⊥

= Z

− (Z

) =

= Z

− (Z

) =

= B

In a recursive approach, the dimension of the ma-

trices involved must not depend on the number

of measurements (N = j − 2i + 1) . Since Q

∈

(m+l)i×(N−2i+1)

, that could be a problem in the

SVD step. However,

) (B

)

= B

(29)

and, if B

= USV

) (B

)

= US

= B

(30)

This means that B

has the same left singular vec-

tors and the same singular values as Z

⊥

. Since

we only need U

and S

to determine a basis for

the column space of Γ

, we only have to deal with

∈ R

li×(m+l)i

, in order to estimate Γ

ICINCO 2004 - SIGNAL PROCESSING, SYSTEMS MODELING AND CONTROL

178

4.2 The system matrices

To estimate matrices A, B, C and D, one must solve

two least squares problems: (4) and (7). In the last

one, both matrices K

and N

have constant di-

mensions, but in the ﬁrst one that does not happen.

Therefore, we will implement (4) recursively, using

the LQ decomposition, to avoid this dimension prob-

lem, since the least squares solution of

i−1

i+1

(i)

ACK

(31)

with

ACK

A K

C K

is the same as the least squares solution of

(N)

= Q

(N)

i−1

i+1

(i)

= (32)

(N)

ACK

To update R

(N)

and B

(N)

we need the last columns

of Y

(i)

, U

, Z

and Z

i+1

. The last column of Z

computed directly from (27) and the last column of

i+1

is computed from

i+1

= Y

f−

H+U−

= (33)

= [

U−

H+

]

(i+1)

where

f−

= L

(i+1)





(i+1)

L3(i+1)





(34)

−

= Y

f−

(i+1)

(35)

H+

= Y

f−

(i+1)

(36)

However a problem may here arise, since the rows-

paces considered at the k-th iteration are not quite the

same as in the (k+1)-th iteration. In order to deal with

this we have implemented two options: recursive least

squares with a forgetting factor or the sliding window

technique.

4.3 The modiﬁed Householder

algorithm

Since our versions of subspace algorithms are se-

quences of least squares problems implemented

through the QR decomposition (or the similar LQ),

we modiﬁed the Householder algorithm for matrix

factorization updating, to updates both B and R, at

once. In fact, we apply the same Householder re-

ﬂection to change both R

(N)

and B

(N)

. As R

(N)

already a upper triangular matrix, the QR decompo-

sition of R

(N+1)

is faster to compute, because now

there is only one line to set to zero. Therefore, the

Householder Algorithm can be adapted in order to

modify only two rows of the matrix in each iteration k

(the last row and the k-th one), since the others remain

the same. Another version of this modiﬁed algorithm

considers the update of the decomposition matrices

when Nr > 1 rows are added.

4.4 Recursive R-MOESP-2QR

algorithm

1. LQ decomposition and the orthogonal projections

(a) Update the last columns of the block Hankel ma-

trices U and Y

U(N +1)

and φ

Y (N+1)

(b) Update the decompositions (18) and (34). Al-

though L

(N+1)

and B

(N+1)

have one more line,

after the LQ decomposition, this last column of

(N+1)

has only zero elements and it can be re-

moved. So, the matrices L

(N+1)

and B

(N+1)

still have the dimensions of L

(N)

and B

(N)

and z

i+1

, the last columns of Z

and Z

i+1

, as in (27) and (33).

2. The column space of Γ

and Γ

i−1

(a) Determine

= U

1/2

, where the system or-

der n is assumed to be known (3) and

H(N+1)

= USV

= U

(37)

(b) Compute

i−1

l(i−1)

3. Matrices A, C and K

(a) Update the QR decomposition matrices R

AC(N )

and B

AC(N )

N+1

= Q

AC(N +1)

= Q

AC(N +1)

AC(N )

N+1

where

N+1

i(N+1)

[

0 I

] φ

U(N +1)

N+1

i−1

i+1(N+1)

[

0 I

] φ

Y (N+1)

(b) Solve, in the least squares sense, the problem

AC(N +1)

= R

AC(N +1)

, where

AC(N +1)

(N+1)

A(N+1)

C(N +1)

RECURSIVE MOESP TYPE SUBSPACE IDENTIFICATION ALGORITHM

179

4. Matrices B and D

(a) Determine K

C(N +1)

, from K

(N+1)

, and matrix

, as in (6)

(b) Solve the least squares problem (7) recursively

K(k−1)

N(k−1)

¸·

(k)

4.5 Recursive PO-MOESP-2QR

algorithm

1. LQ decomposition and the orthogonal projections

(a) Update the last columns of the block Hankel ma-

trices U and Y

U(N +1)

and φ

Y (N+1)

(b) Update the decompositions (18) and (34).

2. The column space of Γ

and Γ

i−1

(a) Determine

= U

1/2

, where the system or-

der n is assumed to be known (3) and

H(N+1)

= USV

= U

(b) Compute

i−1

l(i−1)

0 I

l(i−1)

3. Matrices A and C

(a) Determine

A =

i−1

(b) Determine

C = [

]

4. Matrices B and D (”minimum simulation error

method”)

(a) Update the QR decomposition

BD(N)

¸·

where B

BD(N)

∈ R

m(n+l)×1

, R

BD(N)

∈

m(N+l)×m(N+l)

and

k−1

t=1

⊗ (CA

k−1−t

) ∈ R

l×mn

= u

⊗ I

∈ R

l×lm

= y

∈ R

l×1

= vec(B) ∈ R

mn×1

= vec(D) ∈ R

ml×1

(b) Determine vec(B) and vec(D), by solving,

through a back substitution process,

BD(N+1)

= R

BD(N+1)

vec(D), respectively.

5 RECURSIVE 1QR APPROACH

In this approach, Z

, Z

i+1

and Z

⊥

are obtained

from a unique LQ decomposition:

(i)

0 0 0

0 0

Y 1

Y 2

Y 3













(38)

In fact,

= U

)

−1

= Q

(39)

(40)

and, since the column space of

f−

spanned by the columns of

(i)

−

(41)

Therefore, if we deﬁne

= [

] = (42)

(i)

−

then

= [

]

(43)

i+1

= [

]

(44)

and Z

⊥

= B

5.1 Recursive R-MOESP-1QR

algorithm

1. LQ decomposition and the orthogonal projections

(a) Update the last columns of the block Hankel ma-

trices U and Y

U(N +1)

and φ

Y (N+1)

(b) Update the decomposition (38)

and z

i+1

, the last columns of Z

and Z

i+1

, as in (43) and (44).

2. Determine

= U

1/2

, where

2(N+1)

= USV

= U

(45)

and the system order n is assumed to be known (3).

Compute

i−1

l(i−1)

3. and 4. Implemented as in the algorithm described

in section (4.4)

ICINCO 2004 - SIGNAL PROCESSING, SYSTEMS MODELING AND CONTROL

180

5.2 Recursive PO-MOESP-1QR

algorithm

1. LQ decomposition and the orthogonal projections

(a) Update the last columns of the block Hankel ma-

trices U and Y

U(N +1)

and φ

Y (N+1)

(b) Update the decomposition (38)

2. The column space of Γ

, Γ

and Γ

i−1

(a) Determine

= U

1/2

, where the system or-

der n is assumed to be known (3) and

2(N+1)

= USV

= U

(46)

(b) Compute

i−1

l(i−1)

0 I

l(i−1)

3. and 4. Implemented as in the algorithm described

in section (4.5)

6 NUMERICAL EXAMPLE

It was considered a ”real-life” process (van Overschee

and de Moor, 1996), a laboratory setup acting like

a dryer. The order is four and it can be determined

through the singular values, shown in ﬁgure 1.

Figure 1: Singular values of the sytem, for the ﬁrst 200 mea-

surements. We have assumed that the order of the system is

There is one input only, the voltage over the heat-

ing device (a mesh of resistor wires), and one output,

the air temperature, measured by a thermocouple at

the output – air is fanned through a tube and heat-

eds at the inlet. We have compared only the 1QR

versions, since they produce basically the same es-

timates than the 2QR versions, using i = 15, a forget-

ting factor of 0.995 in the second least squares prob-

lem of the R-MOESP and, for initialization, we have

used the ﬁrst 200 samples of a set of 1000. Figure 2

shows a snapshot of the estimated eigenvalues of A,

with both algorithms and the ofﬂine version of PO-

MOESP (N=200 to N=650). They are all very close

and converging to the true values. In ﬁgure 3, we

show the frequency response of the estimated system,

for the ofﬂine R-MOESP algorithms (left) and the re-

cursive R-MOESP-1QR algorithm (right), for the ﬁrst

350 measurements.

Figure 2: Snapshot (N=200 to N=650) of the trajectory

of the estimated eigenvalues, for the algorithms here pro-

posed (version 1QR) and the ofﬂine version of PO-MOESP

(N=200 to N=650).

Figure 3: Frequency response of the ofﬂine estimated sys-

tem (left) and the recursive estimated system (right).

As to the numerical efﬁciency, the main differ-

ences between the proposed algorithms lie in the ﬁrst

step. In fact, when dealing with LQ decompositions,

a small difference in the size of two matrices has a

big effect on the number of ﬂops required to complete

a LQ decomposition: approximately

+ 2n

ﬂops (Golub and van Loan, 1996), on a n

×(n

+ 1)

matrix (the size of matrices in the online algorithms).

If we compare the less efﬁcient algorithm from Van

Overschee and De Moor (n

= 2 (m + l) i = 180)

with its most efﬁcient version, the R-MOESP-1QR

= (2m + l) i + l = 138), we can see that less

RECURSIVE MOESP TYPE SUBSPACE IDENTIFICATION ALGORITHM

181

than 23% in the number of rows produces less 57.7%

of the number of ﬂops (7.8 · 10

in the Van Over-

schee and De Moor’s algorithm and 3.31 · 10

in the

R-MOESP-1QR version).

7 CONCLUSIONS

In this paper some algorithms for subspace on-line

identiﬁcation have been introduced. They are based

on the PO-MOESP (Verhaegen, 1994) and R-MOESP

(van Overschee and de Moor, 1996) techniques, and

therefore implemented through LQ decompositions.

However, unlike the original algorithms, the proposed

recursive algorithms are based on a least squares in-

terpretation of the orthogonal projections. In fact, Z

i+1

and even Z

⊥

are related to least squares

problems. This allows us to deal with LQ decompo-

sitions of smaller matrices, improving the numerical

efﬁciency of the algorithms without any loss of accu-

racy. We can also use a modiﬁed Householder algo-

rithm, specially developed to improve the efﬁciency

of this LQ-based least squares problems.

Two versions were proposed: the 2QR, where two

LQ decompositions are needed to compute Z

, Z

i+1

and even Z

⊥

, and the 1QR version, where only a

single LQ decomposition is needed. This last version

is, as expected, more efﬁcient (and therefore more in-

teresting) than the ﬁrst one. However, when compared

with the iterative version developed directly from the

ofﬂine algorithm, the 2QR may also be more efﬁcient

than the later. This happens due to the size of the ma-

trices involved. In fact, the LQ decomposition of a

× (n

+ 1) matrix (the size of matrices in the on-

line algorithms), needs approximately

+ 2n

ﬂops (Golub and van Loan, 1996). This means that

two LQ decompositions of lower dimension matrices

may be more efﬁcient than one single LQ decomposi-

tion, of a matrix with bigger n

Further improvements will be made, in order to

consider the update (although partial) of the singu-

lar value decomposition. It is also expected the more

accurate study of the convergence question.

REFERENCES

Chiuso, A. and Picci, G. (2001). Asymptotic variance of

subspace estimates. In Proc. CDC’01 - 40th Con-

ference on Decision and Control, Orlando, Florida,

USA, pages 3910–3915.

Cho, Y., Xu, G., and Kailath, T. (1994). Fast recursive iden-

tiﬁcation of state space models via exploitation of dis-

placement structure. Automatica, 30:45–59.

Delgado, C. and dos Santos, P. L. (2004a). Analysis and

efﬁciency improvement of the estimation of matrices

a and c in some subspace identiﬁcation methods. In

Proc. SICPRO’05 – 4th International Conference on

System Identiﬁcation and Control Problems, Moscow,

Russia (Accepted for Publication).

Delgado, C. and dos Santos, P. L. (2004b). Exploring the

linear relations in the estimation of matrices b and d in

subspace identiﬁcation methods. In Proc. ICINCO’04

– 1st International Conference on Informatics in Con-

trol, Automation and Robotics, Set

ubal, Portugal.

dos Santos, P. L. and de Carvalho, J. L. M. (2003). Improv-

ing the numerical efﬁciency of the b and d estimates

produced by the combined deterministic-stochastic

subspace identiﬁcation algorithms. In Proc. CDC’03

– 42nd Conference on Decision and Control, Maui,

Hawaii, USA, pages 3473–3478.

dos Santos, P. L. and de Carvalho, J. L. M. (2004). A new

insight to the a and c matrices extraction in a MOESP

type subspace identiﬁcation algorithm. In Proc. Con-

trolo’04 – 6th Portuguese Conference on Automatic

Control, Faro, Portugal.

Golub, G. and van Loan, C. (1996). Matrix Computations.

Johns Hopkins University Press, 3rd edition.

Gustafsson, T., Lovera, M., and Verhaegen, M. (1998).

A novel algorithm for recursive instrumental variable

based subspace identiﬁcation. In Proc. CDC’98, Con-

ference on Decision and Control, pages 3920–3925.

Larimore, W. (1990). Canonical variate analysis in identiﬁ-

cation, ﬁltering and adaptive control. In Proc. CDC’90

– 29th Conference on Decision and Control, pages

596–604.

Ljung, L. (1983). Theory and Practice of Recursive Identi-

ﬁcation. The MIT Press.

Ljung, L. (1987). System Identiﬁcation: Theory for the

User. Prentice-Hall.

Moonen, M., de Moor, B., Vandenberghe, L., and Vande-

walle, J. (1989). On- and off-line identiﬁcation of lin-

ear state space models. International Journal of Con-

trol, 49:219–232.

Oku, H. and Kimura, H. (2000). Recursive algorithms of

the state-space subspace system identiﬁcation. Auto-

matica.

van Overschee, P. and de Moor, B. (1996). Subspace Iden-

tiﬁcation for linear systems. Kluwer Academic Pub-

lishers.

Verhaegen, M. (1994). Identiﬁcation of the deterministic

part of MIMO state space models given in innovations

form from input-output data. Automatica, 30:61–74.

Viberg, M., Ottersten, B., Wahlberg, B., and Ljung, L.

(1991). Performance of subspace based state-space

systems identiﬁcation methods. In Proc. IFAC’91 –

12th IFAC World Congress, Sydney, Australia.

ICINCO 2004 - SIGNAL PROCESSING, SYSTEMS MODELING AND CONTROL

182