Frequency Domain System Identication with IV Based Subspace Algorithm

Frequency Domain System Identication with IV Based Subspace Algorithm

Tomas McKelvey

Dept. of Electrical Engineering, Linkoping University S-581 83 Linkoping, Sweden,

Phone: +46 13 282461, Fax: +46 13 282622 Email: tomas@isy.liu.se.

February 27, 1995 Submitted to 34th CDC

Report LiTH-ISY-R-1775, ISSN 1400-3902

Abstract

In this paper we discuss how the time domain subspace based iden- tication algorithms can be modied in order to be applicable when the primary measurements are given as samples of the Fourier transform of the input and output signals. Particularly we study the PI-MOESP algo- rithm 19] in a frequency domain framework. We show that this method is consistent if a certain rank constraint is satised and the frequency domain noise is zero mean and have bounded covariance. An example is presented which illuminates the theoretical discussion.

Keywords: Identication Subspace Method Stochastic Analysis Singular Value Decomposition.

1 Introduction

Methods which identify state-space models by means of geometrical properties of the input and output sequences are commonly known as subspace methods and have received much attention in the literature. The early subspace iden- tication methods 3, 15, 18] focused on the deterministic systems with errors represented at the outputs. Such models are also known as Output-Error Models (OEM). By extending these methods, consistent algorithms have been obtained when the errors are described by colored noise 17, 19, 20]. In spite of the in- herent complexity of these algorithms, many successful applications have been reported. One of the advantages with the methods is the absence of a paramet- ric iterative optimization step. In classical prediction error minimization 11],

ThisworkwassupportedinpartbytheSwedishResearchCouncilforEngineeringSciences

such a step is necessary for most model structures. An excellent overview of time domain subspace methods is given in 21].

In this paper we consider the case when data is given in the frequency do- main, i.e. when samples of the Fourier transform of the input and output signals are the primary measurements. In a number of applications it is com- mon to t models in the frequency domain 16, 12]. Subspace based algorithms formulated in the frequency domain has appeared recently. A frequency domain version of 10] (which is closely related to the basic projection algorithm 3, 2]) has been described in 9] and analyzed with respect to consistency in 13, 14].

A dierent route is followed in in 13] where an algorithm which is based on the inverse discrete Fourier transform and a variation of the realization algorithm

7] is presented. This algorithm is consistent for a large class of systems but re- quires samples of the frequency response of the system at equidistant frequencies covering the whole frequency axis (0

).

The topic of this chapter is to demonstrate that the recent so called sub- space based identication algorithms, by simple modications, can be used to identify systems using frequency domain data. In contrast to the algorithm in

13] wherein only frequency response data was considered here we will assume the Fourier transform of the input and output signals to be given at arbitrary frequencies. The algorithm which we will present is based on the time domain version called PI-MOESP 19]. We will show that the algorithm is strongly consistent for mild assumptions on the noise.

2 Problem Description

Given samples of the discrete time Fourier transform of the input and output signals of a dynamical system we seek an algorithm which identify a state-space model of nite order.

2.1 System Assumptions

Let us consider stable time-invariant discrete time systems of nite order n . One form of describing such a system is by the state-space equations pair

x ( k + 1) = A

x ( k ) + B

u ( k )

y ( k ) = C

x ( k ) + D

u ( k ) + n ( k )  (1) where u ( t )

^m , y ( t )

^p and x ( t )

ⁿ . n ( t )

^p is the noise term which we assume is independent of the input sequence u ( t ). Here the time index k denotes normalized time. Hence y(k) denotes the sample of the output signal y(t) at time instant t = kT where T denotes the sample time. We also assume that the state-space realization (1) is minimal which implies both observability and controllability 6]. A system with this type of noise model is commonly known as output-error models 11]. Note that all such pairs (1) describing the same input/output behavior of the system are equivalent under a non-singular similarity transformation T

^{n n} 6], i.e the matrices ( T

ATT

BCTD ) will be an equivalent state-space realization.

The frequency response of the system (1) is dened as

G ( e ^j! ) = C ( e ^j! I

A )

B + D: (2)

3 Subspace Identication

3.1 The Basic Relations

By introducing

y q ( k ) =

2

6

6

6

4

y ( k ) y ( k + 1) y ( k + ... q

1)

3

7

7

7

5

(3)

and conformally u q ( k ) and n q ( k ), the extended observability matrix

O

q =

2

6

6

6

4

CA C CA ... ^q

3

7

7

7

5

(4)

and the lower triangular Toeplitz matrix

; q =

2

6

6

6

4

D 0 ::: 0

CB D ::: 0

... ... ... ...

CA ^q

B CA ^q

B ::: D

3

7

7

7

5

(5)

we will by recursive use of (1) obtain 2] the relation

y q ( k ) =

q x ( k ) + ; q u q ( k ) + n q ( k ) : (6) The extended observability matrix has a rank equal to the system order n if q

n since the system is minimal.

The discrete time Fourier transform

of a sequence f ( k ) is dened as

F

f ( k ) = F ( ! ) =

k

f ( k ) e

^j!k (7) where j = sqrt

1. From the denition (7) it follows immediately that the discrete Fourier transform of a time shifted sequence satises

F

f ( k + n ) = e ^j!n F ( w ) : (8) Let Y ( ! ) =

y ( k ), U ( ! ) =

u ( k ) and N ( ! ) =

n ( k ). If we now apply the Fourier transform

on both sides of (6) we obtain

W ( ! )

Y ( ! ) =

q X ( ! ) + ; q W ( ! )

U ( ! ) + W ( ! )

N ( ! ) (9) where

denote the Kronecker product 5] and

W ( ! ) =

1 e ^j! e ^j

^!

e ^j!

^q

^T : (10)

Notice that (9) is the frequency domain version of (6).

Assume we have samples of the Fourier transform of the input and output at M frequencies ! = !

:::! M . By collecting these samples in matrices

Y qM =

W q ( !

)

Y ( !

)  W q ( !

)

Y ( !

) 

 W q ( ! M )

Y ( ! M )

^{qp M} 

U qM =

W q ( !

)

U ( !

)  W q ( !

)

U ( !

) 

 W q ( !

)

U ( ! M )

(11)

^{qm M} 

N qM =

W q ( !

)

N ( !

)  W q ( !

)

N ( !

) 

 W q ( ! M )

N ( ! M ) (12)

^{qp M} 

X M =

X ( !

)  X ( !

) 

 X ( ! _M ) 

^{n M}  and using (9) we arrive at the matrix equation

Y qM =

q X M + ; q U qM + N qM : (13) Notice that an alternative expression for U qM can be formulated as

U qM =

W ( !

) W ( !

) ::: W ( ! M )

I m diag( U ( !

) U ( !

) :::U ( ! M )) : In the following discussion we will suppress most size indicating subscripts in (14) order to simplify the notation.

3.2 Identication

The identication scheme we employ to nd an state-space model ( ^ A B ^{^} C ^{^} D ^{^} ) is based on a two step procedure. First the relation (13) is used to consistently determine a matrix ^

q with a range space equal to the extended observability matrix

q . From ^

q it is straight forward to nd ^ A and ^ C as is well known from the time domain subspace methods 21]. In the second step ^ B and ^ D are determined by solving the minimization problem

B ^ D ^{^} = argmin _BD

^M

k

j

Y ( ! k )

( D

C ^{^} ( e ^i!

I

A ^{^} )

B ) U ( ! k )

(15) which has an analytical solution since the transfer function G is a linear function of both B and D assuming ^ A and ^ C are xed.

3.3 The Basic Projection Method

First consider the noise free case and we restate the basic projection method

3, 2] in the frequency domain. In (13) the term ; q U can be removed by the use of 

which is the orthogonal projection onto the null-space of U ^,

^{= I}

U ^{H (} UU ^{H )}

U ⁽¹⁶⁾

here U ^H denotes the complex conjugate and transpose of the matrix U . The inverse in 16 will exist if the system is suciently excited by the input and we return later to this issue. Since

U

_U

= 0

the eect of the input will be removed and we obtain

Y

_U

=

q X

_U

: (17)

Provided

rank( X

U

) = n (18)

Y

_U

^and

q will span the same column space. In the frequency domain for- mulation a small complication occurs at this stage. The state-space matrices ( ABCD ) and thus also the extended observability matrices are usually real matrices but Y

_U

is a complex matrix. The real space can however be re- covered by using both the real part and the imaginary part in a singular value decomposition 4]

;

Re( Y

_U

) Im( Y

_U

)

=

U s U o

 s 0 0  o

V _Ts V _To

(19) where U s

qp n contains the n principal singular vectors and the diagonal matrix  s the corresponding singular values. In the noise free case  o = 0 and there will exist a nonsingular matrix T

^{n n} such that

O

q = U s T:

This shows that U s is an extended observability matrix ^

of the original system in some realization. From U s we can proceed to calculate A and C as

A ^ = ( J

U s )

J

U s (20)

C ^ = J

U s (21)

where J i are the selection matrices dened by

J

=

I

_q

_p 0

_q

_{p p}

 J

=

0

_q

_{p p} I

_q

_p

(22)

J

=

I p 0 _p

_q

_p

(23)

and I i denotes the i

i identity matrix, 0 i j denotes the i

j zero matrix and X

= ( X ^T X )

X ^T denotes the Moore-Penrose pseudo-inverse of the full rank matrix X . With the knowledge of ^ A and ^ C , ^ B and ^ D are easily determined from (15).

3.3.1 Eective Implementation

A most eective way of forming the matrix Y

is by use of the QR factor- ization of the matrix 18]

U Y

=

R

⁰ R

R

Q ^H

Q ^H

: (24)

Straight forward calculations reveal that

Y

⁼ R

Q ^H

and the column space of R

is equal to the column space of Y

^{and it}

suces to use R

in the SVD (19).

3.3.2 Relation With Other Projection Methods

The frequency domain method described in 8, 9] is closely related to the ba- sic projection method presented above. Extend the Fourier transform samples U ( ! k ) and Y ( ! k ) with their corresponding negative frequency value

U (

! _k ) = U ( ! _k )

 Y (

! _k ) = Y ( ! _k )

and form U ^and Y including both negative and positive frequencies and form the projection matrix 

_U

. By then determine the observability range space from the matrix

Y

U

Y ^{H (25)}

we end up with the method described in 9]. Comparing with (17) we conclude that the method of 8] and the basic projection algorithm outlined above will give identical estimates of the observability range space.

3.3.3 Consistency Issues

As we have seen the basic projection algorithm will estimate a state-space model which is similar to the original realization in the noise free case. If we now let the noise term N ( ! ) be a zero mean complex random variable the issue of consistency becomes important. Does the estimate converge to the true system as M , the number of data, tends to innity? Consistency of the basic projection algorithm have been investigated in 13]. The basic projection algorithm is consistent if the frequency data is given at equidistant frequencies covering the entire unit circle and the noise N ( ! k ) is zero mean and have equal covariance proportional to the identity matrix for all frequencies. The uniform covariance requirement for all frequencies and the need for an equidistant frequency grid limits the practical use of the basic projection algorithm.

3.4 Instrumental Variable Techniques

The strict noise properties required in order to guarantee consistency is a severe drawback for the basic projection method. The origin of the problem stems from the fact that the noise inuence does not disappear from the estimate but is required to converge to an identity matrix. What we would like is to nd some instruments which are uncorrelated with the noise but preserves the rank condition (18). The time domain instrumental variable technique in 19, 20] will here be adopted to yield a frequency domain instrumental variable method.

Partition Y , U and N as

Y =

Y p

Y f

U =

U p

U f

N =

N p

N f

The partition is done such that each sub-matrix will have  block rows and consequently 2  = q block rows. The size requirement of this partition is that

 > n where n is the system order. It is straight forward to show that

Y f =

 X f + ;  U f + N f (26) where X f is given by

X f =

e ^j!

X ( !

) e ^j!

X ( !

) ::: e ^j!

X ( ! M )

:

We now have the possibility to remove the future inputs by a projection and then use the past inputs in order to remove the noise 19]. Using (26) we yield

Y

U Hp =

 X f 

U Hp + N f 

U Hp (27) where 

f

denotes the orthogonal projection onto the null-space of U f and is given by

= I

U Hf ( U f U Hf )

U f : (28) If we assume N ( ! _k ) to be zero mean independent random variables with uni- formly bounded second moments

E N ( ! k ) N ( ! k ) ^H = R ( ! k )

R

! k

the following relation follows from a standard limit result 1, Theorem 5.1.2]

M lim

1

M N f 

U Hp = 0  w.p. 1 : If certain conditions are fullled then

rank( 1 M X f 

U Hp ) = n: (29) In the time domain setting this corresponds to certain assumptions of the exci- tation signal 19]. If the rank condition (29) is fullled then the n principal left singular vectors of

Y f 

U Hp

will constitute a strongly consistent estimate of the range space of the extended observability matrix (4).

3.4.1 Implementation

Just as for the basic projection algorithm an ecient implementation involves an QR factorization of the data matrices. By following 19] we form the QR factorization

0

@

U f

U p

Y f

1

A

=

0

@

R

0 0

R

R

0

R

R

R

1

A 0

@

Q ^H

Q ^H

Q ^H

1

A

(30)

By using (30) it is straight forward to show that

Y

U Hp = R

R ^H

and we will use R

in a SVD to estimate the range space of the observability matrix since R

is of full rank whenever U has full rank which we assume. The observability range space is thus extracted as U s

p n from

;

Re( R

) Im( R

)

=

U s U o

 s 0 0  o

V _s ^H V _o ^H

: (31)

Notice that the orthogonal matrix Q in (30) is not needed in the estimation and

the QR factorization constitutes a considerable data reduction since the size of

R

^{p m} is independent of the number of data samples M . As before we use ^ U s as the estimate of the extended observability matrix and determine A and ^ C according to (20) and (21) while ^ B and ^ D are determined from (15).

By using the ^ A and ^ C from the consistent estimates of the observability range space the solution of ^ B and ^ D from (15) is a linear function of the output Fourier transforms Y ( ! k ) and hence also in the noise. By similar arguments as before this shows that ^ B and ^ D will be also will be consistently estimated.

We summarize this discussion in the form of an identication algorithm.

Algorithm 1

1 Form the matrices Y ^{(11) and} U (12) and partition them as

Y =

Y p

Y f

 U =

U p

U f

such that each sub-matrix have  > n block rows.

2 Calculate the QR factorization

0

@

U f

U p

Y f

1

A

=

0

@

R

0 0

R

R

0

R

R

R

1

A 0

@

Q ^H

Q ^H

Q ^H

1

A

3. Calculate the SVD of R

;

Re( R

^{) Im(} R

⁾

= U _s  s V _Ts + U _o  o V _To

were U _s contain the left singular vectors of the n dominating singular values.

4. Determine ^ A and ^ C : A ^ = ( J

U s )

J

U s

C ^ = J

U s

5. Solve the least-squares problem for ^ B and ^ D : B ^ D ^{^} = argmin _BD

^M

k

j

Y ( ! _k )

( D

C ^{^} ( e ^i!

I

A ^{^} )

B ) U ( ! _k )

We also summarize the theoretical discussion in the following theorem.

Theorem 1 Assume that the following conditions are satised:

(i) The frequency data is generated by a stable linear system G ( z ) of order n . (ii) rank( U ) = 2 m

(iii) rank(

X f 

U Hp ) = n

(iv) The noise N ( ! k ) are zero mean independent random variables with bounded covariances

EN ( ! k ) N ( ! k ) ^H = R k

R <



! k

Let ^ G ( z ) be the resulting transfer function when applying Algorithm 1. Then

M lim

sup _z

=1

k

G ( z )

G ^ ( z )

F = 0  w.p 1

Remark 1 From the construction of U (14) we notice that in the single input case m = 1 the rank condition (ii) is equivalent to require that at least 2  samples of U ( ! k ) are non-zero.

4 Illustrating Example

This section describes an identication example based on simulated data. From the results of the example we will clearly see the limits of the basic projection algorithm when faced with data which do not comply with the assumptions needed for consistence. On the other hand the instrumental variable algorithm we will experience to perform as predicted by the consistency result of Theo- rem 1.

4.1 Experimental Setup

Let the true system G ( z ) be a fourth order system with an output error noise model H ( z ). In the frequency domain we thus assume

Y ( ! ) = G ( e ^j! ) U ( ! ) + H ( e ^j! ) E ( ! )

where Y ( ! ), U ( ! ) and E ( ! ) are the Fourier transform of the time domain quantities outputs y ( t ), inputs u ( t ) and innovations e ( t ). The system G ( z ) is given by

G ( z ) = C ( zI

A )

B + D with

A =

0

B

B

@

0

:

:

;0

:

:

0 0 ;0

:

:

0 0 ;6

:

:

1

C

C

A

 B =

0

B

B

@ 0

:

0

:

0

:

0

:

1

C

C

A

C =

^:

^:

^:

^:

 D =

:

:

The noise transfer function is of second order and is given by H ( z ) = C n + ( zI

A n )

B n + D n

with A n =

0 : 6296 0 : 0741

;

7 : 4074 0 : 4815

 B n =

0 : 0370 0 : 7407

C n =

1 : 6300 0 : 0740

 D n = 0 : 2000 :

The Fourier transform of the noise E ( ! ) is modeled as a complex Gaussian dis- tributed random variable with unit variance and is assumed to be independent over dierent frequencies. In the output error formalism we obtain the output error as

N ( ! ) = H ( e ^j! ) E ( ! )

Average

G

G ^{^}

M Proj. Alg. IV Alg.

30 1.7886 1.0425 50 1.3829 0.7737 100 1.2867 0.5078 200 1.2638 0.3751 400 1.2378 0.2550

Table 1: Monte Carlos simulations comparing the basic projection algorithm and the IV algorithm. The estimation error decreases for an increasing amount of identication data which is predicted from Theorem 1. The projection al- gorithms fails to capture the true system which show that the assumption of evenly spaced frequencies and equal covariances are essentially necessary for the projection algorithm to be consistent.

which thus is a complex Gaussian random variable with frequency dependent variance equal to

H ( e ^j! )

. The Fourier transform of the input signal is dened to be U ( ! ) = 1 

! , i.e. all frequencies are equally excited.

To examine the consistency properties of the basic projection algorithm and the IV algorithm we perform Monte Carlo simulations estimating the system given samples of U ( ! ) and Y ( ! ) using dierent noise realizations of E ( ! ) and an increasing number of samples of the transforms. The frequency grid will be logarithmically spaced between !

= 0 : 3 and ! M = . Data lengths of 30, 50, 100, 200 and 400 frequency samples will be used. For each data length 100 dierent noise realizations are generated and both algorithms estimate 100 models. To assess the quality of the resulting model the innity norm of the estimation error

k

G ( z )

G ^ ( z )

is determined for each estimated model and averaged over the 100 estimated models.

4.2 Estimation Results

As expected from the analysis the quality of the estimates from the instrumen- tal variable algorithm (IV) improves as the number of samples of the Fourier transform increases. In Table 1 the averaged maximum identication error is presented. The results clearly indicate that the basic projection algorithm is not consistent for these data. We have in this example violated the requirement of equally spaced frequencies and equal noise covariances required for consistency of the basic projection algorithm 13] and by judging from the example these requirements seems to be essentially necessary for consistency.

5 Conclusions

In this paper we have shown that the PI-MOESP instrumental variable algo-

rithm can be modied in order to be applicable when the primary measurements

are given as samples of the Fourier transform of the input and output signals.

Error Alg. 10.2 Error Alg. 10.1 G(z) H(z)

0 0.5 1 1.5 2 2.5 3 3.5

10⁻⁴ 10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹

Frequency (rad/s)

Magnitude

Figure 1: Result from Monte Carlo simulations using data length M = 400. The true transfer function G ( z ) is depicted as \+" and the noise transfer function H ( z ) is shown as the dotted line. The absolute value of the mean transfer function errors calculated over 100 estimated models are shown as a solid line for the IV method and as a dashed line for the basic projection method.

We have shown that the method is consistent if a certain rank constraint is satis-

ed and the frequency domain noise is zero mean and have bounded covariance.

An example is presented which illuminate the theoretical discussion.

References

1] K. L. Chung. A Course in Probability Theory. Academic Press, San Diego, CA, 1974.

2] B. De Moor. Mathematical Concepts and Techniques for Modeling of Static and Dynamic Systems. PhD thesis, Katholieke Universiteit Leuven, Kard.

Mercierlaan 94, 3030 Leuven (Heverlee), Belgium, June 1988.

3] B. De Moor and J. Vandewalle. A geometrical strategy for the identi- cation of state space models of linear multivariable systems with singular value decomposition. In Proc. of the 3rd International Symposium on Ap- plications of Multivariable System Techniques, April 13-15, Plymouth, UK, 1987.

4] R. A. Horn and C. R. Johnsson. Matrix Analysis. Cambridge University

Press, Cambridge, NY, 1985.

5] R. A. Horn and C. R. Johnsson. Topics In Matrix Analysis. Cambridge University Press, Cambridge, NY, 1991.

6] T. Kailath. Linear Systems. Prentice-Hall, Englewood Clis, New Jersey, 1980.

7] S. Y. Kung. A new identication and model reduction algorithm via singu- lar value decomposition. In Proc. of 12th Asilomar Conference on Circuits, Systems and Computers, Pacic Grove, CA, pages 705{714, 1978.

8] K. Liu, R. N. Jacques, and D. W. Miller. Frequency domain structural system identication by observability range space extraction. Technical report, Space Engineering Research Center, MIT Cambridge, MA 02139, September 1993. Submitted for review for 1994 ACC and ASME Journal of Dynamic Systems, Measurement and Control.

9] K. Liu, R. N. Jacques, and D. W. Miller. Frequency domain structural system identication by observability range space extraction. In Proc. of the American Control Conference, Baltimore, Maryland, volume 1, pages 107{111, June 1994.

10] K. Liu and R. E. Skelton. Q-markov covariance equivalent realization and its application to exible structure identication. AIAA Journal of Guid- ance, Control and Dynamics, 16(2):308{319, 1993.

11] L. Ljung. System Identication: Theory for the User. Prentice-Hall, En- glewood Clis, New Jersey, 1987.

12] L. Ljung. Some results on identifying linear systems using frequency do- main data. In Proc. 32nd IEEE Conference on Decision and Control, San Antonio, Texas, pages 3534{3538, December 1993.

13] T. McKelvey and H. Ak#cay. An ecient frequency domain state-space identication algorithm: Robustness and stochastic analysis. In Proc. 33rd IEEE Conference on Decision and Control, Lake Buena Vista, Florida, pages 3348{3353, December 1994.

14] T. McKelvey, H. Ak#cay, and L. Ljung. Ecient construction of transfer functions from frequency response data. Technical report, Report LiTH- ISY-R-1609, Dep. of EE, Link$oping University, S-581 83 Link$oping, Swe- den, 1994.

15] M. Moonen, B. De Moor, L. Vandenberghe, and J. Vandewalle. On- and o- line identication of linear state-space models. Int. J. Control, 49(1):219{

232, 1989.

16] J. Schoukens and R. Pintelon. Identication of Linear Systems: a Practical Guideline to Accurate Modeling. Pergamon Press, London, UK, 1991.

17] P. Van Overschee and B. De Moor. N4SID: Subspace algorithms for the

identication of combined deterministic-stochastic systems. Automatica,

30(1):75{93, 1994.

18] M. Verhaegen. A novel non-iterative MIMO state space model identication technique. In Proc. 9th IFAC/IFORS Symp. on Identication and System parameter estimation, Budapest, Hungary, pages 1453{1458, July 1991.

19] M. Verhaegen. Subspace model identication, Part III: Analysis of the ordinary output-error state space model identication algorithm. Int. J.

Control, 58:555{586, 1993.

20] M. Verhaegen. Identication of the deterministic part of MIMO state space models given in innovations form from input-output data. Automatica, 30(1):61{74, 1994.

21] M. Viberg. Subspace methods in system identication. In 10th IFAC

Symposium on System Identication, volume 1, pages 1{12. IFAC, 1994.