Order Selection of Models for Time Domain Identification of Transfer Functions

CO(>\ ri g III ~ . !\ li lll.tlit)1l

CD

IF. \(

I ~ I ;-.i .-l .

I dt'II1 ili( ,lIi olJ ,(lid

Ytlrk .

s\ . . 1t'1I I

P ; II -. IIll t' l tT

L· K . I q;-.i .-)

ORDER SELECTION OF MODELS FOR TIME DOMAIN IDENTIFICATION OF TRANSFER FUNCTIONS

",

Yuan Z-d* and Shen L-j** "' n e/HI , 11111 '111 of .\fll lfll' lIIlllin. LIl"I Chi/Ill .\'omllli [ '/I;"('/'" ily. Sfllll/gfllli.

r hl' I'mpil''' UI'/)//b/i( of Chil/Il :;. :/) ('/)(11"111/(' 111 "/ .\/ otlif'lIwtirol .' ilotistin. FfI., 1 ChiN" .\'unll"/ l' lIi1. '(' ,-.\ ily. SIIfIII,eJlfli.

Thl' Pm/JIt"" RI'/mlt/i( of Chil/Il

Abstract.

Based on Ljung and Yuan's paper(198?), the problem of select-

ing model orders for indentification of transfer functions is considered. The rate of increase of deN) gives a conservative estimate of the model order.

For given accuracy of the transfer function estimate, we

can

determine the order d(N) and data length N by calculating the residual terms.

Order selection using hypothesis testing is also discussed.

KeyWords.

1• In

Transfer functions; stochactic system; Identification.

INTRODU CTION ~jung

The model set is taken to be d : yet) = Gd(e , q-1)u(t) +v(t)(1.3) J'ii d

and Yuan's paper (1985), the

transfer function is parametrized as a

where d

Gd(~' q-1) =L.

black box and no given order is chosen a priori.

This means that the model

d t)

= (gl' .... 'gd)

( 1 .4)

'1'

If we introduce

number of observed data tends to infi-

fdT(t) = (u(t-1), ... ,u(t-d»)

From this point of view, some

the model can be rewritten yet) = fdT(t)e d + vet)

interesting properties of the transfer function estimate have been derived under certain assumptions.

-k

and

order may increase to infinity when the

nity.

gk • q

k=1

( 1.6)

the well known least square estimate

Following

R (N)-1 1

the notations used in the Ljung and

N

L.

N t=1

d

fd(t)y(t)

where

Yuan's (1985) paper. ConSider a stochastic system S

: yet)

where

= Go (q-1)u(t)

+ vet)

( 1 .8)

The estimate of the transfer function is A ieA) .. (1 .9) Gl'J (e ) = Wd * here

(1.1)

-

-k ( 1 .2) q g~ k=1 -1 and q is the delay operator, vet) a

Go ( q -1)

(1. 7)

L

~

Wd =

(-ie,.)

e

, e

-i2 e,.)

, ••• ,e

-idc.>

I

( 1 .10)

List below 1jung and Yuan's results(1985):

zero mean stochastic process.

As

I

I ~)

I

N ___

and deN) _

0...

I I ~I:!

YU
(a )

Consistency sup

-1T~<.U~iT

I'2N(e i

'-»

!

~ (ei "') _ ..~ ( i N .t'."N e

LV

AS N ---..

C>
The conditions of Ljung and Yuan's

As ymptotic variance expression N

Tllli IL,'l't; 0F INCFU:;ASE OF d(N)

2.

_Go( c ic..:> j! _ O w.p.1.

(b)

L~ i

J2 ~

theorems on the model order can be summarized

as the following (Ljung and

Yuan 1985): as k --+ 00 (2.1 )

1)1 :

lJ2 :

Asymptotic normality

(c)

d(k) ---0-"'"

E k=1

lJ3 :

as k --..

l d(k 2) /k J 2<

+

00

(2.2 )

()o

d(k) non-decreas1ng ; AsN (0,

1) I.V !

d(k)/k

1»c.:»)

non-inc reasinp, . 3 [d (k)/k JP<'+"for some P>O

D4 :

Two questions on the model order could be asked: (1)

Nhy should we assume that d(N)

(2)

---. be as N--+ 00 ?

What is the optimal order

(2.4 ) It.,..,1s well known that if

L. N=1

.1 < + 00 NO<

integers and

~>

1 then

where N tak es positive .

,

0(

real numbers.

,

With

selection as a function of the data

this result we can easily prove the

length

following Lemma.

N'!

The second problem can be partly solved

Lemma

2.1

o

by finding out the rate of increase of d(N) as N--.. oo.

Ljung and Yuan (1985)

have given a rough formula for selecting the model order d(N):

Proof.

shall

give

d(N) =,O(NO< ),

r

,some constant,

O.

01..>

Examine condition D1 and D3.

LO ~ (N'" )

d(N)

a

N

conservative estimate of d(N) that may

J

)'20(N2 0(

2

N

N

that is

be more useful.

d 2 (N)

For given accuracy of the transfer function estimate, to find the suitable d(N) and N is another problem discussed in sect ion 3.

< 1/4,

D2 is obvious because

2

In section 2, we

< 0<

then D1, D2 and D3 hold.

d(N) = O(log N),

i.e., lim(d(N)/logN) = constant.

I f d(N) =)"O(N'" ),

0

2

• 0 (

)

N1 -20(

N Due to

0(

< 1/4, 1-2<>« 0, then

2 d (N)/N--+ 0

The model order

as N

--~

selection using hypothsis testing is given in section 4.

Some numerical

=0 0

i l lustration and conclusions are given in section 5 and 6 resp ectively.

and

2-40(

>

1 (--) N2 - 4 o(

1,

hence C>o

L N=1

<

00

)

Id e ntifi ca tion of' Transkr FUllctioll s

We can also select d(N) non-decreasing 2 and d (N)jN non-increasing, for example

¥.

d ( N) =

<

< 0(

NO( , 0

U4 :

11 93

C ~

cp u (

as N __

Using this lemma we can obtain a series

T,r N( T) __ • 0

T=-2d (N)

Jd (N)

The proof ls complete.

~ "() > 0

~N )

U5: _1_

1 /4 •

tU )

.0<0

hold, then

of modified theorems that correspond to Ljung and Yuan's theorems ln their paper (1985) • Theorem 2.1

Consider the transfer i~

A

function estimate ~(e


) obtained from

(1.6) - (1.10).

Proof.

Assume that

Yuan (1985) and Lemma 2.1.

n(N)

(11)

, 0 <0«1/4 (2.5)

OO(No/)

U1:

U2:

J

I u( t) l< C, \I t

and

fN( w, d( N»

11minf N

-.00

Theorem 2.4.

m-dependent sequence E(v(t»4 ~ C

(2.6) :?;

'S >

Assume that the conditions

of Theorem 2.3 hold and{v(t)~iS a

the input sequence {u( t) \

ls independent of the nolse sequence { v( t)

Use Theorem 4.1 of Ljung and

0

(2.7)

where

q,v (~ )

d

fN(W, d(N»

t

=.L. r (T)e-i~T T=-d N (2.8) u(t-T)u(t)

(2.9)

N T=1

asN( 0, Proof.

cPU<

c.) )

Use Theorem 5.1 of Ljung and

Yuan (1985) and Lemma 2.1. The rate of increase d( N) = )" O(N 0(

(i11)

V1

Ev(t) = 0 /Ev(t)v(S)

(2.10)

I<

C1 A t-StO<~(1 (2 .11 )

then

I ~N( e iW )

sup

- EGN( e i

c.,) )

I

-"* 0

w.p.1 • Proof .

gives an estimate of the model order which guarantees that the theoremS hold. For example,

N

= 1000, 4IN ~ •. = 5.62,

hence d = 6.

Use Theorem 2.1 of Ljung and DETERlUNI NG d(N) AFD N

Yuan (1985) and Lemma 2.1. Theorem 2.2 00

L.

FuR GIVEN ACCURACY

If the conditions of Theo-

ConSider the accuracy of the transfer

rem 2.1 hold and d ( N)

I g~ \--...

k=d+1

0

as d (N) --.....

function estimate. Let

then

B(d(N» = EGN(eic.,) w.p.1. Proof.

as

N-_

Use Theorem 3.3 of Ljung and

Theorem 2.3

( i) t (11), U1, (i11) and

N»

--~fu( <..l)

(GN(eit.» + (i>G,_ ( e i <.>

Assume that the condition

fN( "", d(

- Go(eic..:l) (3.1 )

In view of GN( e iw ) _ Go(ei w )

0..

Yuan (1985) and Lemma 2.1.

U3:

)

l"

_ EGr (eic.,)) iw » ) _ (G ( e

+

0

we define the mean squa re error as N-~~

1'; (0

, Ii • w )

I I ~J4

i'uan Z-cl and She n L-j

with the asymptotic expression for the variance of GN(ei~), we have M( d ,N,

Co)

q,)w)

d(N)

_.!.-

)

L

I

gO,)

k=d+1 + B 2 (d)

k

Jd .. c

vhere

q:u(,,»

N

C>o

~ ~ (C

11

R- 1 (N)11

<\

,Iu(t)l
d

we obtain The MSE i.e. M( d ,N, c..;»

depends on a

IB(d)\~\W. E (e~

e~)\+I~«(..)\

-

number of design variables, such as input spectrum, data length, model structure and model order.

For given

MSE (accuracy) and given bias (can be estimated), ve try to determine the suitable d(N} and N.

Revrite the true

y(t) =

g~

L

q-k

(d

*

2 C

~+

u(t) + v(t)

k=1

Theorem 3.1

f:>

for

I f Theorem 2.3 holds and

0

IB(d)12~(d.c2S+

Introduce

....

k=d+1

<

~t


N

L

(3.4)

k=d+1 Summarizing these calculations, we obtain

d(N) _ _ _4> _ v(w) _

~(t)

I g~ \

L

1)

the folloving theorem.

system (1.1) and (1.2) c>c>

(>0

~

1)(

E

k=d+1

19~\)2
g~ u( t-k) ~

where

+

~

= 1, then

MSE(d, N) Proof.

<[

From (3.3), (3.5) and (3.6),

we obtain (3.7). then Go(e iW ) = W"

«(..) )8~

From (1.7), ve have

+

~ «0.')

For given MSE( i.e .given

~

> 0), we can use

(3.6) to determine d(N) and use (3.5) to determine N.

The numerical illustration

will be given in section 5. In (3.3), the terms are non-negative, and us i ng ( 2 .5), we can obtain upper bound of d(N) +~(

t»

Since sup -Ti~c..>~;r

l'u (w

) (,5.8)

(,.8) also gives a rough es timation of d(N) for given N, NSr.; and B(d).

II ~5

Idelltificatioll of Tralls icr FUllctioll s

4.

ORDER SELECTIUN USING

For simplicity, we use the

HYPOTHESIS TESTING

criterion to replace the above mentioned.

In Theorem 3.1, it is difficult to esti-

~ ollowin~

If

(4.8 )

mate NSE'and B(d) when the true transfer function is unknown. 't:e

11 0W

pothesis testing method.

develop a hyCollect two

groups of data (from two independent ex-

we reject the null hypothesis, we

~ust

increase the model order d 1 • If sup IU(d1'd 2 ,w)l
(4.9)

-lf~t..:>~Tf

From the

periments) of data length N.

we accept the null hypothesis and take

first group, we obtain

~ '" -iku:> ( 4 .1 ) =1~_e k=1 r, From the second ~roup of data, we obtain '" ic..:l G (d , e ) N 1

"" GN( d 2 , e i~ )

"" -ikw = $2 L- gk e

(4.2 ) k=l According to Theorem 2.4 (provided the

5.

NUl
Consider the true system q

y( t)

S1:

-1

_ -1 u( t) + e( t)

conditions of Theorem 2.4 hold for both

1-0.2:>q that is

groups of data) ,

LGr-/d1 ' AsN(O,

yet) = E(1/4)k-\-l
'" (d ' e i~ ) - t;G N 1 cv )

d1
N4>u C ..,j

[

~ G (d , N 2

e

ic..:l

AsN(O, d 1 N

.A

) - t;G ( d , e N

iu:>

2

4>v
output respectively.

) 1E.

a

sequence of independent rectangle (4.4)


f e( tHiS

where d2 > d 1 • The statistical null hypothesis can be written as

distributed random variables over the interval

[-1/2, 1/2 ] •

The input u(t) is given by u(t)

= H(q-1)*

w (t)

where 1

H(q-1)

k

(?Q

k

=2:.(0.3) q1-0 .3q -1 k=O { w( t) }is a sequence of independent rectangle distributed random variables over the interval [-1/2. 1/2J. Under the null hypothesis Ho' U(d 1 , d2'~ ) has an asymtotic normal distribution

1'lOreover, {e( t)

and

{ 'W ( t) J are

independent. Obviously, we can obtain

N( 0 ,1) •

lu(t)j

If the confidence level is given, we can obtain the critical value

1

f~

means we must continue to increase the model order d 1 • I f IU(d 1 ,d 2 , ...,)1 < f,,("'), we accept the null hypothesis, that means we can take deN) = d 1 •

....

f:o(0.3)k

~ ~0(0.3)k

so that

p fIU(d1, d2 , )\ ~ fo< 1= 1 - 0( (4.7) I f the actual value Iu( d 1 ,d , w ) j ) f 0( 2 then we reject the null hypothesis, which

~

and

f

~r..,(t-k)1 ; 0.5 = 5/7

~ g~ (1/4)k-1 k=d+1 . k=d+1 Consider the system S2: u(t) = H(q-1)

it

we have

Tu (w ) and (k( CA> )

1/241T

wet)

11%

i'uan Z-d and Shen

L~j

(N = 1000, d

hence

6) and G (e i "') o


MIP

)

1.0') - 0.6 L:os (I.)

;,:4'IT

-

Let Rd(N) be the d)( d covariance matrix

- - - - - ---

--

of {u(t)} (0.5) I j-il 1-( 0 .5)2

12 PHASE

Simple calculations show that Rd = 12 •

1+(0.5)

-U.5

2

l

dxd

w.

o

_

i<.J

) - - --

•

Steyart, 1973), the eigenvalue of

R~1 satisfies ~ 12 .. 0.5

The conservative estimate of deN) gives

or

a rough estimate of the model order.

I A - 12 .. (1 + 0.3 2 ,) I ~ 12 .. 2 .. 0.3

By calculating the mean square error

Hence

and the bias, ye can verify the model

/I.~ 12 .. (1+0.3)2 = 1.69 .. 12 = ~ 'l.'he norm

11

CONCLUSIONS

6•

I

I'>.. - 12

order and data length, which guarantees

R~1" is defined as

that the theorems hold.

"R~11\= max

{Ai i=1 ..• d -1 Aiare the eigenvalues of Rd

1

(1000~:: 5.62), from (5.10)

Taking d = 6 we obtain that

2 (d .. C

b + 1 ) 2 ( L. I g~ \

k=d+1 2 25/49" 12" 1.69 + 1) • )( (1/3 2 ) (1 /4 2 (6-1 ») = 0 .00042

and 6

Asymptotic properties of black-box identification of transfer functions. IEEE

Transaction of Automatic

Control, June 1985. (To appear.) (1.09-0 .6 Cos(w»

Astrom, K. J. (1970).

1000 6

.. 1.09

Introduction

to stochastic control theory,

0.01014

Academic Press.

1000 Vue to (3.8) we obtain

Stewart, G. W. (1973).

hSE(d,N) ~ 0.01014 + 0.00042

Introduction

to matrix computations, Academic

0.0106 "0.02

'I'herefore, for

HE~'ERENCES

Ljung, L. and Z. D. Yuan. (1985)

)2

- (6

~

The Bode plots of GN( e

and G (e i <.»

According to Gerschgorin Theorem (G.

A

Fig. 1.

Press. 1000

=0 .0 2 we have N

and d = 6. A

Figure 1 Shows the Bode plots of GN(e

ic.l

)