An Ecient Frequency Domain State-Space Identication Algorithm:

An Ecient Frequency Domain State-Space Identication Algorithm:

Robustness and Stochastic Analysis

Tomas McKelvey and Huseyin Akcay Department of Electrical Engineering

Linkoping University S-581 83 Linkoping, Sweden e-mail:

Fax: +46 13 282622 February 23, 1994

33rd CDC

Abstract

In this paper we present a novel algorithm for identifying linear time-invariant discrete time state-space models from frequency response data. The algorithm is non-iterative and exactly recovers a true system of order

, if

+ 2 noise-free uniformly spaced frequency response measurements are given. Analysis show that if the measurements are perturbed with errors upper bounded by

the identication error will be upper bounded by

and hence the algorithm is robust. An asymptotic stochastic analysis show, under weak assumptions, that the algorithm is consistent if the measurements are contaminated with noise. In a companion paper the algorithm is applied to real data with promising results.

Keywords: system identication, state-space methods, frequency response.

1 Introduction

It is often desired to t parametric models to experimental data in order to gain information about a system in a condensed form. Models of state-space type are favorable since modern control design methods, and simulation tools are based on state-space models. Recently new sub-space algorithms have emerged which identies state-space models from time-domain input/output data, see 13, 15, 10]. These algorithms are related to the realization algorithms by Ho and Kalman,

5] and Kung, 8] or the ERA-algorithm 6], which all identify/realize a state-space model given

the impulse response (the Markov parameters) of the system. However there are cases in which

frequency response data are available rather than time-domain data. With the advent of modern

spectral analyzers and high performing data-acquisition equipment, which produce high quality

frequency response data, it is of practical interest to nd techniques which identify high-quality

state-space models from frequency data.

Classically, a parametric model has been tted to the frequency response by some curve tting technique. An early result for SISO models in the Laplace s domain is 9]. This approach gives an implicit high frequency weighting which makes the algorithm biased. In 14] an iterative scheme (SK-iterations) were suggested to remove this deciency. This approach has however some severe convergence problems as pointed out in 17].

To directly obtain state-space models, the Inverse Discrete Fourier Transform (IDFT) have been used in order to acquire estimates of the impulse response coecients and then proceed by using a realization algorithm, see e.g. 7]. The disadvantage with this approach is the distortion of the estimates introduced by the IDFT (time-aliasing), since, in reality, only a nite number of frequency responses can be used.

In order to circumvent this problem Bayard, see 1], introduced an algorithm (SSFD) which, as a rst step, ts an overparametrized rational model to the data by iteratively minimizing (by SK-iterations) a 2-norm error criterion and then determine the impulse response from the rational model. A state-space model is then estimated by the ERA-algorithm, see 6].

This paper introduce and analyze a novel frequency domain state-space identication algorithm using the frequency response of a system. The major features of this new approach are: A system of order n is exactly recovered using n + 2 frequency response samples if the data is noise-free.

The method is non-iterative. Analysis demonstrates that the algorithm is robust to norm bounded perturbations. Furthermore if the frequency response measurements are corrupted with stochastic noise the algorithm is consistent (the true model is recovered as the number of measurements tend to innity). The algorithm involves 3 key steps. First an IDFT is performed to obtain (distorted) impulse response estimates. Secondly the estimates are used in a realization step, by a singular value decomposition (SVD), to obtain the A and C matrices. Even though the impulse response estimates are distorted these estimates of A and C are correct (in the noise-free case). In the third step B and D are estimated by an ordinary least-squares method. In a companion paper the algorithm is applied to real data, originating from a exible mechanical structure, with promising results.

2 Problem formulation

We will assume that the true system G is a stable multivariable linear time-invariant discrete time system with input/output properties characterized by the impulse response coecients g k through the equation

y ( t ) =

k

g k u ( t

k ) (1)

where y ( t )

IR ^p , u ( t )

IR ^m and g k

IR ^p

^m . If the system is of nite order n it can be described by a state-space model

x ( t + 1) = Ax ( t ) + Bu ( t )

y ( t ) = Cx ( t ) + Du ( t )  (2) where y ( t )

IR ^p u ( t )

IR ^m  and x ( t )

IR ⁿ . The state-space model (2) is a special case of (1)

with g k =

(

D k = 0

CA ^k

B k > 0 : (3)

The frequency response of (1) is

G ( e ^j! ) =

k

g k e

^j!k (4)

which for the state-space model (2) can be written as

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

A )

B + D: (5) The problem formulation is then Given a nite number M + 1, possibly noisy, samples

G k = G ( e ^j!

) + e k  ! k

0  ] (6) of the frequency response of the system nd a nite dimensional state-space system (2) of order n , denoted by ^ G , such that the true system and the identied model are \close", where the closeness is quantied by the following distance between the true and estimated transfer functions

jj

G

G ^{^}

= sup

_!   ( G ( e ^j! )

G ^ ( e ^j! )) : (7) Here   ( A ) denote the largest singular value of the p

m matrix A . Furthermore, denote by  i ( A ) the i th singular value, where the ordering is given by 

( A )



( A )

:::

 r ( A )

0 with r = min( pm ).

3 The algorithm

If the impulse response coecients (3) are given, well-known realization algorithms can be used to obtain state-space realizations, see e.g. 5] and 8]. The algorithm to be presented is closely related to these results but does not require the true impulse response coecients to be known.

Assume that frequency response data G ( e ^j!

) on a set of uniformly spaced frequencies, ! k =

k M  k = 0 ::: M are given. Since G is a real transfer function, frequency response data on 0  ] can be extended to 0  2  ] as follows

G ( e ^j

^k=M

) = G

( e ^j

^M

^k

^=M )  k = 1 :::M

1 (8) where (

)

denotes complex conjugate. Construct the impulse response block Hankel matrix

H ^ qr =

2

6

6

6

6

4

g ^

g ^

::: g ^ r

g ^

g ^

::: g ^ r

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

g ^ q g ^ q

::: g ^ q

r

3

7

7

7

7

5

(9) with number of block rows q > n and block columns r

n and where the rst q + r

1 impulse response coecients of G are estimated from the 2 M -point Inverse Discrete Fourier Transform (IDFT) according to

g ^ i = 1 2 M

2

M

k

G ( e ^j

^k=

^M ) e ^j

^ik=

^M  i = 0 :::q + r

1 : (10)

The size of ^ H qr is limited by q + r

M .

Compute the singular value decomposition of ^ H qr

H ^ qr = ^ U

U ^{^}

]

"

^

0 0 ^

#"

V ^

^T V ^

^T

#

(11) where ^

contains the n principal singular values.

The system matrices are then estimated as

A ^ =

I

_q

_p 0

_q

_p

_p

U ^

^

⁼

0

_q

_p

_p I

_q

_p

U ^{^}

^

⁼

(12) C ^ =

I p 0 _p

_q

U ^

^

⁼

 (13) B ^ D ^{^} = arg min

B

D

M

X

i

G ( e ^j

^i=

^M

)

G ^ ( e ^j

^i=

^M

)

(14) where

G ^ ( z ) = ^ D + ^ C ( zI

A ^{^} )

B ^ (15) and I i denotes the i

i identity matrix and 0 i

j denotes the i

j zero matrix and A

= ( A ^T A )

A ^T denotes the Moore-Penrose pseudo-inverse of the full-rank matrix A . Notice that B and D appear linear in (15) and thus (14) can be solved by a linear least-squares solution.

From (10), notice that

M lim

g ^ k =

^

0

G ( e ^j

^ ) e ^j

^k d = g k  k = 0 :::q + r

1  (16) where g k is the k th impulse response coecient of G . Let H qr and  i  U i  V i  i = 1  2 denote the limits of the matrices in (11{15) as M tends to innity. The resulting algorithm is then known as Kung's algorithm 8], in which ^ B is calculated from

B ^ = 

⁼

V

^T

"

I m

0 m

r

#

: (17)

The realization given by (11-15) in Kung's algorithm 8] is balanced in the sense that the q -block row observability matrix

O

q =

2

6

6

6

6

4

CA C CA ... ^q

3

7

7

7

7

5

(18)

and the r -block column controllability matrix

C

r =

B AB ::: A ^r

B

(19)

satises

O

Tq

q =

r

Tr = 

(20)

A true balanced realization will thus be obtained only if both q and r tend to innity. In this limiting case 

equals the Hankel singular values of the system G , see 3]. For nite q and r the Hankel singular values of G are underestimated by 

. Although the singular values 

do not play any role (except for a selection of a base for the state-space variables) in the construction of a state-space model for G they will be essential in the selection of a model order in the presence of unmodelled dynamics and noise. (Notice that A and C can be calculated from (11-13) letting

= I n ). As q and r increase, the singular values of H qr also increase in the absence of noise.

However, small singular values of H qr and corresponding left and right singular vectors will be less reliable in the presence of noise. Therefore q and r should be chosen suciently large to obtain 

as large as possible and the model order should be kept as low as possible. We will formalize these observations in a succeeding section.

The spectral radius of a square n

n matrix is dened as

( A ) = max

:

( A )

(21)

where ( A ) are the eigenvalues of A .

The usefulness of the above presented algorithm in the case of nite M follows from the following key theorem.

Theorem 3.1 Let G be an n th order system represented by (2). Assume that = ( A ) < 1 . Then n +2 noiseless equidistant frequency response measurements of G on 0  ] are sucient to identify G by the above algorithm (8{15) and

^

(1 +

^M )

: (22)

Proof : Since G is a stable transfer function, it can be represented by the following Taylor series G ( z ) = D + C ( zI

A )

B = D +

k

CA ^k

Bz

^k =

k

g k z

^k  <

z

: (23) Notice that ^ g k can be written as

^ g k =

i

g ( k + 2 iM ) = CA ^k



1

X

i

A

^iM

!

B = CA ^k

( I

A

^M )

B = CA ^k

B ^~ (24) and therefore ^ H qr can be factored as

H ^ qr =

q

~ _r =

q  ~ B A B ::: A ^{~ r}

B ^~ ] : (25) Hence, if r

n and q > n , then ^

= 0 and the column range spaces of H qr and ^ H qr will be equal.

Thus ^ A C ^{^} G ^{^} are related to AC and G by a square nonsingular matrix S as

A ^ = S

AS  ^ C = CS  ^ G = ^ D + C ( zI

A )

S B: ^{^} (26) Since (14) has a unique solution, we get ^ B = S

B , ^ D = D: Then ( ABCD ) and ( ^ A B ^{^} C ^{^} D ^{^} ) are similar. Letting q = n + 1  r = n M = n + 2, we satisfy the condition q + r < 2 M .

For the second part: notice that ^ H qr H ^ _Tqr

(1+

^M )

H qr H _Tqr is positive semidenite if ~

r

~ _Tr

(1 +

^M )

r

Tr is positive semidenite which is true.

In the algorithm given by (11{14), the singular values of ^ H qr only inuence the basis chosen for the state-space variables and not the transfer function G . Since in a realistic problem, ^

is never zero due to noise and unmodelled dynamics, a clear separation between the singular values of ^

and ^

must be obtained in order to choose a model order unambiguously. Theorem 3.1 gives a bound on the distortion of the smallest singular value due to aliasing caused by the nite 2 M -point inverse DFT.

4 Analysis

In the last section, we demonstrated that an n th order system can be recovered by a simple al- gorithm involving only the IDFT, SVD, and a linear least-squares solution using n + 2 frequency response measurements. The assumption that the system which generated the data is nite di- mensional with a known order is not realistic and we will admit the possibility of some \small unmodelled dynamics" which is not captured by any nite dimensional model set. However, by Theorem 3.1, the model order can be selected quite large and by choosing a suciently large model order, one may hope to reduce unmodelled dynamics. A competing factor to unmodelled dynamics is statistical uncertainty of the model due to the disturbances which grows as model order increases.

In this section, we will demonstrate that the proposed algorithm is robust to unmodelled dynamics and noise. In the following section we analyze the inuence of noisy data and show consistency for the algorithm using mild assumptions on the noise. For complete proofs of the results below we refer to 12].

4.1 Unmodelled Dynamics

Suppose that we are given frequency response data

G k = G ( e ^jk=M ) + e k  k = 0  1 ::: M (27) where G represents the transfer function of a n th order stable, single-input/single-output (SISO) nominal model and e k captures unmodelled dynamics assumed to be uniformly bounded

k

e

: (28)

We present two important perturbation lemmas which are used in the derivation of transfer function error bounds. The rst lemma gives an upper bound on the singular values of H qr whose elements are constructed from e k .

Lemma 4.1 Let G = 0 and let frequency response measurements of G be corrupted as in (28).

Then _M lim

!1

  ( ^ H qr )

1 + log( q + r

1)] : (29) The bound in the lemma is accurate even for a modest number of data. It is stated in the asymptotic form for convenience only. If the Hankel structure of ^ H qr and the fact that ^ H qr is constructed by the IDFT is not exploited only a conservative bound

(( q + r )  ) could be derived.

Let

A

denote the range space of the matrix A . Let G be n th order and

U

denote its

observability range space obtained by the SVD of H qr as in (11) when M tends to innity. We will

bound small perturbations of U

in the Frobenius norm caused by the unmodelled dynamics and

the DFT. Let

U ^

denote perturbed observability range space.

Lemma 4.2 ^Let G be an n th order system (2) with = ( A ) < 1 . Suppose that the frequency response of G is corrupted as in (28). Consider the SVD of H qr obtained as in (11). Assume that

 _n

( H qr ) > 2 (1 +

^M )

2 q

2 (1

^M )

G

+ 

 (30) where = 1 + log( q + r

1) . Then there exists a unique matrix P satisfying

k

P

F

(1 +

^M )

2 q

2

G

(1

^M )

+ 

 _n

( H qr )

(1 +

^M )

2 q 2

G

(1

^M )

+  ]  ⁽³¹⁾ such that

U

+ U

P

=

U ^

for all suciently large M .

We have expressed error bounds and a condition in terms of Mnqr

G

 and in the proceeding discussion, we observed how the choice of M a#ected the singular values and the ob- servability range space. To proceed our discussion , we make the following simplications

  ( ^ H qr ) =   ( H qr )

log( q + r )  ( when G = 0) (32)

 _n

( H qr ) > 4

2 q log( q + r )

G

  (33)

k

P

F

2

2 q log( q + r )

G

 _n

( H qr )

2

2 q log( q + r )

G

: ⁽³⁴⁾ The condition (33) suggests the following bound on

P

F

k

P

F

2

2 q log( q + r ) 

_n ( H qr )

G

: (35) We state the main result of this section in the following theorem

Theorem 4.3 ^Let G be an n th order system (2) with = ( A ) < 1 . Suppose that frequency response of G are corrupted as in (28). Consider the SVD of H qr obtained as in (11). Let G ^ denote the identied model by the Algorithm (8{15). Assume that

 _n

( H qr ) > 2 (1 +

^M )

2 q

2 (1

^M )

G

+ 

 (36) Then the the following transfer function error bound approximately holds

k

G

G ^{^}

16 (1

)

G

n

P

F : (37) Suppose that we model the true system by a nite-dimensional linear-time invariant system with a xed order n . In this case, our lack of information or will to neglect the detailed features of the system such as small time-variations, nonlinearities etc. are reected by the quantity e having a bounded norm. Then letting qr

( n ), Theorem 4.3 shows that

G

G ^{^}

(  ) : Thus the algorithm proposed in the previous section is robust to perturbations bounded in norm.

Now, assume that the order of the nominal model is not predictable a priori , and such that it

is to be inferred as a part of the identication procedure. (In what follows, unmodelled dynamics

is embedded in G ). If certain smoothness conditions are imposed on the system, then it can be

approximated by low order models. Consider, as an example, the class of strictly stable systems.

A system is said to be strictly stable if its impulse response satises

1

X

k

k

g k

<

: (38)

As q and r tend to innity, H qr converges to the Hankel operator for G , see citeGlover:84, and the algorithm (9{15) then gives an n th order truncated balanced realization for G: Thus in the limit, truncation errors for balanced realizations hold for identication errors, i.e.,

k

G

G ^{^}

= 2

k

n

 k ( G )  (39)

where  k ( G )  k = 1  2 ::: denote the Hankel singular values of G .

The error bound in Theorem 4.3 depends on  n ( H qr ) : Assuming G is strictly stable,  n ( H qr ) can be related to the n th Hankel singular value of G as follows

j

 n ( H qr )

 n ( G )

k

min

qr

j

g k

o ( q

) ( q

r ) : (40) which also hold for all singular values  i ( H qr )  i = 1 :::n . Hence for strictly stable systems the singular values of H qr approach the Hankel singular values of G as q

.

4.2 Asymptotic analysis in the stochastic case

In this section we investigate the asymptotic properties of the presented algorithm for the case when the transfer function measurements are obtained from a system (2) of nite order n and each measurement is corrupted by noise. We limit ourself to perform the analysis for the single input multiple output case to keep notation as simple as possible. The same results, however, trivially holds also for the multiple input case. The analysis given below is somewhat related to the time-domain analysis in 16].

Assume that the noisy transfer function measurements are

G ~ k = G ( e ^j!

) + e k (41) where each component of the noise term is assumed to have the form

e k =  k + j k (42)

and  k and  k are assumed to be zero mean independent Gaussian variables with equal variance E k  _Tk = E  k  _Tk =

R k such that

E e k e _Hk = R k (43)

which is the denition of a complex Gaussian random variable, see 4]. In (43) (

) ^H denotes the

complex conjugate and transpose of a complex matrix. We also assume that E e k e _Hl = 0 

l

= k ,

i.e. the noise term for di#erent frequencies are independent. These noise assumptions are rather

weak and, for example, valid if the frequency response is obtained as the empirical transfer function estimate, see 11]. Let  R be the smallest matrix satisfying

R 

R k 

k (44)

where the inequality A

B means that A

B is positive semidenite.

Recall algorithm (8{15). If we construct the Hankel matrix using the noisy samples of the frequency response we obtain the Hankel matrix

H ~ qr = ^ H qr + E qr (45)

where the noise Hankel matrix is built up by the elements

 E qr ] _lm = 1 2 M

2

M

k

e k e

^

^l

^m

^k=

^M : (46) We construct a covariance matrix

~

H

= ~ H qr H ~ _Tqr (47)

which preserves the range space of ~ H qr and we are interested in investigating the properties of this matrix as M

. The left singular vectors of ~

converge w.p. 1 to

H

= lim _M

!1

~

H

(48)

and by using the dominated convergence theorem 2], we obtain

H

= lim _M

!1

E

^~ : (49)

Using the independence of the noise and the transfer function we nd

E

^~ = ^ H qr H ^ _Tqr + E

E qr E _Tqr

: (50) The  lm ] block of the second term is

h

E E qr E _Hqr

_lm = E ₄ M ¹

r

X

k

0

@ 2

M



e  e ^j

^

^l

^k

^=

^M

1

A



0

@ 2

M



e _H e

^j

^

^k

^m

^=

^M

1

A

= r

4 M

2

M



e ^j

^

^l

^m

^=

^M E

E  E _H

(51) which is norm bounded as

jj

h

E E qr E _Hqr

_lm

r

2 M

R 

(52)

Since ^ H qr

H qr as M

we arrive at

H

= lim _M

!1

~

H

= H qr H _Tqr (53)

if the noise covariance is bounded. The left singular vectors of

then span the same space as the observability range space of the true system and hence ^ A and ^ C will be consistent estimates (up to a similarity transformation). Since ^ B D ^{^} are estimated by a linear least-squares solution these estimates will also be consistent. We have thus established the following result.

Theorem 4.4 The estimates ( ^ A B ^{^} C ^{^} D ^{^} ) from algorithm (8{15) converges w.p. 1 to a true system description as M

if the measurement noise has bounded covariance and r and q are kept nite.

This result is quite powerful since consistency is preserved even if noise is correlated between di#erent outputs and/or if the variance is unequal over di#erent frequencies. The latter case arises when the measurement noise is colored in the time domain 11]

If we impose a further restriction on the noise

R k = R

k (54)

the noise term (51) will simplify to

h

E E qr E _Hqr

_lm =

(

r

4

M

R l = m

0  l

= m ⁽⁵⁵⁾

This is interesting since the factor

_M ^r

now allows us also to let q and r , as well as M tend to innity while still preserving consistency if we let

rM lim

r

4 M

^{= 0} : (56)

If the true system G is innite dimensional this allows us to use the known bound (39) for approx- imation by a balanced realization since in the limit, rqM

, we obtain an n th order balanced realization from algorithm (8{15). We formalize this into a theorem.

Theorem 4.5 Assume the noise satises (54). The algorithm (8{15) then, w.p. 1, yields an n th order balanced realization ( ^ A B ^{^} C ^{^} D ^{^} ) as qrM

in a way such that lim rM

r

4

M

= 0 . Furthermore the identication error is bounded by

jj

G

G ^{^}

2

k

n

 k ( G ) (57)

where  k ( G ) are the Hankel singular values of the true system G .

4.3 Discussion

For practical use of the algorithm it is important to know how the choice of qr and M inuence the result. The results above give us some guidelines. If the noise characteristics are unknown we have to rely on Theorem 4.4 which requires M to be larger than qr to average out the noise. This is however conicting with the wish to obtain 

as large as possible since 

is enhanced if qr increase. However 

grows faster than 

, the singular values of the noise, and for a suciently large qr and large enough M a clear separation will occur between 

and 

and the estimated model of order n will be reliable.

If we can assume the noise properties in Theorem 4.5 the choice is clear. Use as large Hankel

matrix as possible to obtain a near-balanced realization and thus get a good rational approximation

to a possible innite dimensional system.

5 Conclusions

We have in this paper presented a non-iterative frequency domain state-space identication algo- rithm. If the frequency data is noise free and is generated by a n th order system we show that only n + 2 equidistant frequency samples are required to exactly recover the true system. If the measurements are contaminated with unmodelled dynamics bounded by  we show that the result- ing identication error is upper bounded by  and hence the algorithm is robust. An asymptotic stochastic analysis shows that the algorithm is consistent if each measurement is perturbed by an independent stochastic noise term.

References

1] D. S. Bayard. An algorithm for state-space frequency domain identication without windowing distortions. In Proc. of the 31st IEEE Conference on Decision and Control, Tucson, Arizona , pages 1707{1712, December 1992.

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

3] K. Glover. All optimal hankel-norm approximations of linear multivariable systems and their L

-error bounds. International Journal of Control , 39(6):1115{1193, 1984.

4] N.R. Goodman. Statistical analysis based on a certain multivariate complex gaussian distri- bution (an introduction). Annals Math. Stat. , Vol. 34:152{176, 1963.

5] B. L. Ho and R. E. Kalman. E#ective construction of linear state-variable models from in- put/output functions. Regelungstechnik , 14(12):545{592, 1966.

6] J. N. Juang and R. S. Pappa. An eigensystem realization algorithm for modal parameter identication and model reduction. J. of Guidence, Control and Dynamics , 8(5):620{627, 1985.

7] J. N. Juang and H. Suzuki. An eigensystem realization algorithm in frequency domain for modal parameter identication. Journal of Vibration, Acoustics, Stress, and Reliability in Design , 110:24{29, January 1988.

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

9] E. C. Levy. Complex curve tting. IRE trans. on automatic control , AC-4:37{44, May 1959.

10] K. Liu and R. E. Skelton. Q-markov covariance equivalent realization and its application to

exible structure identication. AIAA Journal of Guidance, Control and Dynamics , 16(2):308{

319, 1993.

11] L. Ljung. System Identication: Theory for the User . Prentice-Hall, Englewood Cli#s, New

Jersey, 1987.

12] T. McKelvey, H. Ak$cay, and L. Ljung. Ecient construction of transfer functions from fre- quency response data. Technical report, Report LiTH-ISY-xxxx, Dep. of EE, Link%oping Uni- versity, S-581 83 Link%oping, Sweden, 1994.

13] M. Moonen, B. De Moor, L. Vandenberghe, and J. Vandewalle. On- and o#-line identication of linear state-space models. International Journal of Control , 49(1):219{232, 1989.

14] C. K. Sanathanan and J. Koerner. Transfer function synthesis as a ratio of two complex polynomials. IEEE Trans. Automatic Control , 8:56{58, January 1963.

15] M. Verhagen. 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.

16] M. Viberg, B. Ottersten, B. Wahlberg, and L. Ljung. A statistical perspective on state-space modeling using subspace methods. In Proc. 30th IEEE Conference on Decision and Control, Brighton, England , pages 1337{1342, 1991.

17] A. H. Whiteld. Asymptotic behaviour of transfer function synthesis methods. Int. Journal

of Control , 45(3):1083{1092, 1987.