The Use of Realization Theory in the Robust Identification of Multivariable Systems

THE USE OF REALIZATION THEORY IN THE ROBUST IDENTIFICATION OF MUL TIV ARIABLE SYSTEMS G. A. van Zee and O. H. Bosgra Department of Mechanical Engineering, Delft University of Technology, Mekelweg 2, 2628 CD Delft, The Netherlands

Abstract. The ro.bustness o.f predictio.n erro.r identificatio.n metho.ds for appro.ximate mo.dels is discussed and is sho.wn to' depend en the availability o.f initial estimates. The paper provides a so.lutio.n to. this pro.blem by pro.po.sing a numerically stable realizatio.n procedure fer sample co.variance functio.ns o.f measured variables, which is based upo.n a singular value deco.mpo.sitio.n o.f an appro.priate Hankel matrix. Fro.m the realizatio.n thus o.btained o.ne finds an appro.ximate state-space mo.del fer the pro.cess and disturbances under co.nsideratio.n via a facto.rizatio.n o.f the asso.ciated spectrum. The realizatio.n pro.blem is solved in a transparent and ro.bust way, sho.wing directly the po.ssible trade-off between mo.del o.rder and perfo.rmance. Furthermore, the po.sitive-real pro.perty o.f the co.variance realizatio.n is investigated and is sho.wn to. define a lo.wer bound fer the co.variance matrix with argument zero.. This pro.vides a test en the co.nditio.ning o.f co.variance realizatio.ns and suggests a method to. improve ill-co.nditio.ned realizatio.ns. The perfo.rmance and ro.bustness o.f the pro.po.sed metho.d is demo.nstrated by sho.wing numerical applicatio.ns. Keywo.rds. Identificatio.n; multivariable system; minimal realizatio.n; system o.rder reductio.n. 1. INTRODUCTION The estimatio.n o.f multivariable dynamic pro.cess mo.dels for co.mplex industrial plants in the chemical, petro.leum and po.wer industry co.nstitutes a co.nsiderable task fro.m the experimental and co.mputatio.nal point o.f view. Due to. the limited number o.f measured variables in co.mpariso.n with the co.mplexity o.f the physical pheno.mena, o.nly appro.ximate mo.dels fer co.mplete plants can be o.btained fro.m experimental data. Furthermo.re, the intuitively appealing o.bjective of parsimo.ny is in co.nflict with the o.bjective of accuracy, fo.rcing the order of a model to be a co.mpromise as has clearly been expressed by Akaike (1973,1974a) in his information criterion. For these reasons, in several recent papers the identification pro.blem has been considered as an approximatio.n problem (Ljung, 1978; Caines, 1978; Anderson, Mo.ore and Hawkes, 1978; Ljung and Van Overbeek, 1978). A key idea is this line of thought is the co.ncept o.f robustness as a desirable pro.perty of an identification method. In the statistical literature (Huber, 1972) an estimation method is called ro.bust if small deviations fro.m the ideal hypo. theses under which it has been derived cause only small erro.rs in the estimate. Fer identification metho.ds the co.ncept of ro.bustness can be defined similarly. Recently, Ljung (1978), Caines (1978) and Anderson, Moore and Hawkes (1978) have sho.wn that the impo.rtant class o.f 477

predictio.n error identification methods, which includes the maximum likeliho.o.d method and equatio.n erro.r metho.ds, has a certain ro.bustness pro.perty. They show that if the true system is not an element of the parametrized model set, then the o.ptimal model, i.e. the model associated with the glo.bal minimum of the criterion function, co.nverges to. an approximate model that is close to the true system as the number of o.bservations tends to infinity. This general convergence result can replace properties as unbiasedness and co.nsistency that are meaningless in the approximation co.ntext. However, the demonstrated robustness of the optimal model is only of practical value if the global minimum of the criterion can be actually co.mputed. Since the predictio.n erro.r formulation leads to a nonlinear optimizatio.n problem that has no clo.sed fo.rm so.lutio.n, o.ne must reso.rt to a numerical optimizatio.n procedure. Even if such a procedure converges to a so.lution at the expense of co.nsiderable computatio.nal costs, o.ne obtains a local minimum. Experiences in practical applicatio.ns (Eklund and Gustavsson, 1973; Gupta and Mehra, 1974) and results of theo.retical analyses (Bohlin, 1971; Soderstrom, 1974; Astrom and Soderstrom, 1974; Soderstrom 1975) have shown the generic existence of mUltiple lo.cal minima in such problems. A similar phenomenon o.f multiplicity o.f local minima has been observed by Aplevich (1973) in the closely related problem of optimal model approximatio.n.

G. A. van Zee and O. H. Bosgra

478

G.A. van Zee and O.H. Bosgra These observations lead to the point of view that the success of a prediction error identification method depends heavily on the availability of a well defined model set comprising a set of initial model orders, model structures and associated sets of initial parametrizations, each initial parametrization being close enough to the global minimum of the criterion function to assure convergence of a numerical optimization procedure to this minimum. In this approach, the initial parametrization must serve the purpose of a elobal search whereas the prediction error identification formulation provides a final refinement of the model. Several authors have stressed the role of an initial estimate (Tse and Weinert, 1975; Akaike, 1976; Caines and Rissanen, 1974) and have suggested the use of realization theory for an efficient computation of initial model estimates directly from the sample covariance function of the measured variables (Faurre and Marmorat, 1969; Kailath and Geesey, 1971; Rissanen and Kailath, 1972; Rissanen, 1971; Akaike, 1974b). However, there are two main difficulties in such an approach. Firstly, realization theory does not provide numerically well-conditioned procedures for noisy or error-corrupted covariance data. Secondly, a covariance realization constitutes a model for input-output relations and disturbances only if it is positive-real. In this paper, these two problems will be analyzed and it will be shown that for both problems robust and numerically wellconditioned approximate solutions can be obtained. In section 2, the general problem will be considered. In section 3, a singular value decomposition is proposed as the basis of a stable approximate realization procedure. The qualitative and quantitative role of the positive-real condition is investigated in section 4, and in section 5 numerical examples will show how these ideas may lead to robust initial estimates for model structures and parameters. 2. POSITIVE-REAL COVARIANCE FUNCTIONS For a stationary zero mean random proce s s {zK} with a normal distribution, the probability density function p(ZN) of an observation sequence ZN={zO zl ... zN-I} satisfies NZ -! T -I N p~ZN)=«27T) IvNI) exp(-! ZN VN Z ) (2.1)

::::{ZNZNT}.[~~ ~r" '~-Il ~_I

RO

(2.2)

. zT } , z ER Z (2.3) ~ K+~ k Thus the covariance function R. completely determines the observation seq~ence ZN in a probabilistic sense. If the normality assumption doe s not hold, the covariance function still gives an important although not R. =E {z,

complete description of the observations. Let a model for the observation sequence {zk} be given by the innovations representation n Z XER , v,z €R xk+I=A x k + K v k

In the case that the system (2.4) is subjected to a known independent input sequence {uk}, E{uk}=O, E{Uk+iU~}=ROi' then (2.4) becomes

~+I=A ~ zk =

K)[~~

+ (B

[~~] = l~ ]xk

+

[~~J

(2.5)

so that without loss of generality (2.4) can be taken as the model. The problem to be considered is: find a representation ( 2 .4) such that the associated covariance function of the output sequence {zk} is a good approximation to a given sample covariance function Ri of the observation sequence ZN This involves two separate steps: approximating Ri by a finite-order_(loworder) realization, and factoring Ri (i.e. performing a deconvolution) such that the model (2.4) is obtained as a factor . These two steps can be executed in both consecutive orders. Rissanen (1973) and Morf (1974) have presented algorithms to obtain a Cholesky factorization of the block-Toeplitz matrix -T -T

~o

V = ~I N

~I" ·RN-I R O

(2.6)

~_I··· · ··RO

in which the successive block rows of a Cholesky factor converge to values that directly determine the impulse response of an innovations representation. Approximation by a finite-order representation yields a model (2.4). The algorithm requires a numerically well-conditioned VN>O, and estimation errors in Ri that yield oscillatory behaviour for increasing i may cause poor convergence of the algorithm. The reverse orde. of these two steps will be cOGs~dered here: first approximate Ri by a smooth low-order model, and then factor the result to obtain (2.4). Consider the covariance function matrix es timates (Jenkins and Watts, 1968): N-i-I N-i-I I \ T -, I \ T (2.7) R·- N L zk . zk ,R · - N_, L zk+ · zk ~ k=O +~ ~ ~ k=O ~ Definition. A (sample) covariance function Ri, i=O,I,2, •.. is called positive real (p.r , ) on [O,N-I] if VN (2.6) s atisfies VN~ O.

Thus p.r. on [O,N-I] implies p.r. on [O ,N'-I], N' SN. In section4 a p.r. condition on [O,oo ] ... ill be shown to be a necessary condition in order that a finite-order realization of a covariance function can be represented by an innovations model ( 2 .4). Thus both by defining a continuation for a

The use of realization theory

479

Realization theory in robust identification (sample) covariance function such that its domain is extended from [O,N-IJ to [O,ooJ, and by approximating a (sample) covariance function on [O,ooJ, the p.r. condition on [0,00] may be violated, and these two steps must be considered to arrive at a model (2.4). Note that the continuation step alone (e.g. Rissanep. and Kailath, 1972) may yield a result that is not p.r. on [0,00]. These observations underline the importance of a proper choice for a covariance estimator: ~L (2.6) generally does not satisfy a p.r. condition, whereas the (biased) estimate Ri (2.6) is p.r. on [O,N-IJ as shown in appendix I. Evidently, the covariance function of a zero-mean stationary random process is p.r. on [0,00] so that the p.r. condition can be viewed as a structural constraint that a symmetric matrix function must satisfy in order to represent a covariance function. 3. APPROXIMATE COVARIANCE REALIZATION It is well known that the innovations model (2.4) yields a covariance function Ri of the form i-I . B ,~=1,2,.. (3.1) Ri = CA for some matrix B of appropriate dimensions. Let a covariance function estimate Ri, i=l, 2, ... with Ri=O, i~N define the block Hankel

matr[~~ ~~.~~~.~~ Hr=

r

2

rel1ated

:3·····~'+1

,HA=

~ ~+I" ·~+M'-I

[~Il

R ~= ~

~

~

ma[~~ice~3~~:~~~~;~

23

~4""~'+2

~+I ~+2"~+M' ~

,HC=[R I R2"'~']

(3.2)

where the integers M,M' define the dimensions. These matrices constitute a polynomial system matrix (Rosenbrock, 1970) ZHl;-HA HBj (3.3) P(z) = -H

r

C

°

which defines the transfer function matrix G(z)=Rlz-I+R2z-2+ .. if the matrices (3.2) are taken sufficiently large (Bosgra and Van der Weiden, 1978). Define u,r,v by the singular value decomposition

Hr=U~ ~lv,

r=diag(ol,o2,,,.,Ok)

(3.4)

where OI~02~ .. ~ok~0 are the singular values of Hr and the square matrices U,V are orthogonal, i.e. U-I=UT, V-I=VT . Associated with (3.3) is a finite-dimensional linear system with an order equal to rank Hr. A lower-order approximation of (3.3) and thus of the associated sequence RI ,R 2 , •.. can be obtained by extracting from the set of equations implied by the system matrix (3.3) Such linear combinations of equations that these in a definite sense contribute most s~gnificantly to the formal rank of Hr. This wlll be done by virtue of the following theorem. Let rl=diag(ol,o2, •. ,on),let n~rank Hr and suppose that a lower-order

approximation of order n is required. Theorem. Amongst all matrices of rank n with dimensions identical to those of Hr' the matrix

Hi

=U[gl

~]v

minimizes I ~r-Hil IF'

where I 1.1 IF is the Frobenius norm. Proof: see Stewart (1973). 0 Thus by expanding the set of equations (3.3) according to a singular value decomposition of Hr and rejecting linear combinations of equations and variables associated with the smallest singular values, we obtain a loworder approximation in the sense made precise in the theorem. Write U=[UI U2], vT=[vt V~J, where UI,Vt contain n columns. preI and postmultlPl y (3.3) by

° °J

[vir~ (3.5) [rl uT 0] Il' Il respectively, then (3.3) is brought to the

°

;~(~)

= [Z:g-A

~]

(3.6)

with P*(z) of dimensions (n+l)x(n+l) and

_ -! T T-!

A-r l UIHAvlr l (3.7) _ -! T _ T-! B-r l UIH B C-HCvlr l A detailed discussion on these results is given by van Zee and Bosgr~ (1979). Note that the final result (3.7) has been proposed along rather intuitive arguments by Zeiger and McEwen (1974), though their results have not been fully appreciated (Jaaksoo and Nurges, 1975; De Jong, 1978). Note further that in obtaining the solution (3.7) along a singular value decomposition of Hr, which is known to give trustworth judgements about the rank of a matrix, the main goal has been to obtain a robust solution. In section 5 further aspects will be discussed. The result can be summarized as follows. The transfer function matrix associated with (3.6), (3.7), -I G*(z)=C(zIn-A) B (3.8) constitutes an approximation to the transfer function matrix G(z)=RIZ-~R2Z-2+ .. defined by (3.3). Thus a finite-order approximation to Ri, i=0,1,2, .• has been determined as R*=R . 12 (3.9) o 0' R~=CAi-IB ~ , ~= , , •• and by its construction, the matrix triple (A,B,C) constitutes a minimal realization. Moreover, asymptotic stability of A may be assumed because an unstable A must be rejected as an apparent failure of the covariance modeling procedure. 4. INVERSE COVARIANCE GENERATION The relationships between a finite-dimensional covariance function (3.9) and the underlying innovations representation (2.4) are well known (Kalman, 1965; Anderson, 1969; Willems, 1971; Faurre, 1970, 1976). Let the covariance function realization R*i, i=0,1,2, .. (3.9) be denoted by the matrix quadruple (A,B,C,~). Theorem. The covariance function realization

G. A. van Zee and O. H. Bosgra

480

G.A. van Zee and O.H. Bosgra

(A.B.C.~) admits an innovations representation (2.4) if and only if Ri. i=O. 1.2 •..• is positive-real on [O.ooJ. Proof: see Faurre (1976). 0 It is well known that this condition can equivalently be expressed as a condition on the existence of a spectral factor of the associated spectral density matrix or as a condition on the existence of a real symmetric solution to the so-called linear matrix inequality. The following corollary provides a computational test*to investigate a matrix quadruple (A.B.C.Ro). Corollary. The covariance function realization (A.B.C.R;) is positive-real on [O.ooJ if the conditions a and b hold: a. IA(z)i has no zeros on Izl=l. where ZI-~ BR*-IBTj A(z)= CTR:-1C Z-?I-A~ and ~=A-BR:-I~t.l)

[

b. R*+C(I-A)-IB+BT(I-AT)-ICT>O o

(4.2)

Proof: follows directly from Molinari (1975). theorem 2. 0 The corollary is a sufficient condition. and does not consider the pathological cases where IA(z) I. which defines the zeros of the associated spectral density matrix. has paired zeros at the unit circle which might eventually lead to a positive-real result. If IAklfO. the zeros of IA(z) I are the eigenvalues of the Hamiltonian matrix (Molinari. 1975) M=

[

~-BR:-IBT(~:)-ICTR:-IC

_(~)-lcrR:-Ic

The computation of these eigenvalues and the verification of (4.2) provides a basis to assess the positive-real property for (A.B.C.~). If (A.B.C.~) is positive-real. then computation of the innovations model (2.4). i.e. of the matrices K.Q follows from the stabilizing solution of an algebraic matrix Riccati equation (Gevers and Kailath. 1973; Faurre. 1970. 1976; Denham. 1975; and others). If M turns out to have one or more separate eigenvalues on the unit circle. then (A.B.C.~) violates the positive-real condition and must be rejected as a model for a covariance function. This may occur. for instance as the result of choosing a very low or a very high order of the approximation (A.B.C.~). As a logical consequence of the approximations involved in section 3. one may search for a ~~~ which closely approximates ~ but which makes (A.B.C.~) positive-real. The existence of such a R~ is guaranteed by (2.6). In appendix 11. a simple optimization procedure is described to compute ~ such that 11~_~112 is minimized and (A.B.C.IQ;) is positive-real. The experiments in section 5 indicate that violation of the positive-real condition may occur whereas minor modifications of ~ are sufficient to obtain a feasible IQ;. It is an interesting open problem to characterize a (possibly non-unique) lower bound for R6 on the basis of a given (A.B.C).

5. EXPERIMENTS In this section the implications of the outlined problems and the applicability of the proposed solutions are explored by means of simulated experiments on a stirred tank with the variations of the two feed flows with different concentrations as inputs and with the variations of the outcoming flow and its concentration as outputs (Fig. I). A discrete-time 4th order state-space model with system- and measurement noise as described in detail by Kwakernaak and Sivan (1972). with independent zero mean normal distributed white noise as inputs. has been used in a discrete-time simulation. Figure 2b shows three elements of the sample covariance function of a sequence of 8192 observations consisting of the augmented inputs and outputs. A comparison with the exact covariance function in Fig. 2a reveals the estimation errors. Figure 3 shows the singular values of a Hankel matrix HE with dimensions MlxM'l (M=M'=32. l=4). which contains the normalized sample covariance function. i.e. correlation function. An important decision in formulating HE is to choose the dimensions M and M'. These should be large in order to reduce the relative contribution of high frequency estimation errors to HE' On the other hand these dimensions should not exceed certain limits because sample covariance functions are known to be erroneous for large argument values. These considerations lead to a compromise for M.M'. For the given values the singular value decomposition required one minute computing time on an IBM 370/158 computer in PL/I. If the selection of the order of a realization is based upon Fig. 3. it seems reasonable to compute realizations with order 2.3 and 4. because the first three singular values represent almost the complete rank of HE' As expected on the basis of Fig. 3. the accuracy of the realization with order 4 is not better than the third order realization. Therefore only the realizations with order 2 and 3 are plotted in Fig. 4. Both realizations give a smooth and stable approximation of the sample covariance function. but the accuracy of the third order realization is significantly better for element Ry2y2 ' Tests on the positive-real property of the realizations show that the realizations of order 3 and 4 are positive-real. For the second order realization. ~ is adjusted by an optimization procedure as described in appendix 11. Using the Powell (1964) method. the optimum is obtained in 20 seconds computation time. The solution is close to ~ with only a significant change (+5%) in (~)y2y2 as shown in Fig. 5 . Figure 6 shows two elements of the impulse response matrices of the exact innovations model. the second order approximation. and the third

R6

The use of realization theory

481

Realization theory in robust identification order approximation. Table I shows the corresponding intensities Q of the innovations processes. The impulse response elements describing the input-output relations are accurate in both approximate models. The gains of the impulse response elements describing the innovations-output relations show differences that correspond to the errors in Q. These errors are mainly due to different relative heights of the peaks in (R&)y?y2 and (Ro)y2y2, compared to the exact value 1n Ry2y2 (Fig. 5). This part of the approximate models is therefore less accurate than the deterministic part. CONCLUSIONS A solution has been proposed for two interrelated problems that inhibit a straightforward application of realization theory to sample covariance functions. Firstly a globally approximating realization method has been proposed which uses a numerically well-conditioned singular value decomposition. Secondly the positive-real condition has been shown to play a crucial role in the selection of a feasible order of a covariance realization, such that the existence of an innovations representation is guaranteed. It is concluded from simulated experiments that smooth and stable positive-real covariance realizations are easily obtained.Moreove~ if the positivereal condition is violated, e.g. due to a low order of the ap~roximation, then a minor adjustment of RC is sufficient to obtain a positive-real solution. The accuracy of the computed innovations representations seems good enough for the purpose of initial estimates in a prediction error identification method, although further experiments are needed to support this statement. The proposed procedure as a whole is robust and transparent as it provides insight to the user in the conditioning of each step. ACKNOWLEDGEMENT The authors wish to express their gratitude to J. Barendregt, who provided much of the experimental material described in section 5. REFERENCES Akaike, H. (1973). Information theory and an extension of the maximum likelihood principle. In B.N. Petrov, F. Csaki (Eds.), 2nd Int. Symp. on Information Theory, Akademiai Kiad6, Budapest, pp. 267-281. Akaike, H. (1974a). A new look at the statistical model identification. IEEE Trans. Autom. Control, 19, pp. 716-723. Akaike, H. (1974b). StochastTc theory of minimal realization. IEEE Trans. Autom. Control, 19, pp. 667-674. Akaike, H. (1976). Canonical correlation analysis of time series and the use of an

information criterion. In R.K. Mehra and D.G. Lainiotis (Eds.), System Identification: Advance~Case Studies. Academic Press, New York. pp. 27-96. Anderson, B.D.O. (1969). The inverse problem of stationary covariance generation. Journ. Statist. Physics, I, pp. 133-147. Anderson, B.D.O., J.B. Moore and R.M. Hawkes (1978). Model approximations via prediction error identification. Automatica, 14, pp. 615-622. Aplevich, J.D. (1973). On the multiplicity of locally-optimal approximations of linear systems. In Advances in Cybernetics and Systems Research, vol. I. Transcripta, London. pp. 13 17. Astrom, K.J. and T. Soderstrom (1974). Uniqueness of the maximum likelihood estimates of the parameters of an ARMA model. IEEE Trans. Autom. Control, 19, pp. 769-773. -Bohlin, T. (1971). On the problem of ambiguities in maximum likelihood identification. Automatica, 7, pp. 199210. Bosgra, O.H. and A.J.J. van der Weiden (1978). Input-output in.ariants for linear multivariable systems. Report 147-N, Delft Univ. Technol., Lab. Meas. Control, Delft, The Netherlands; submitted to IEEE Trans. Autom. Control. Caines, P.E. and J. Rissanen (1974). Maximum likelihood estimation of parameters in multivariable Gaussian stochastic processes. IEEE Trans. Inform. Theory, 20, pp. 102-104. Caines, P.E. (1978). Stationary linear and nonlinear system identification and predictor set completeness. IEEE Trans. Autom. Control, 23, pp. 583-594. De Jong, L.S. (1975)~Numerical aspects of realization algorithms in linear system theory. Doctoral dissertation, Eindhoven University of Technology, Eindhoven, The Netherlands. Denham, M.J. (1975). On the factorization of discrete-time rational spectral density matrices. IEEE Trans. Autom. Control, 20, pp. 535-537. Eklund, K. and I. Gustavsson (1973). Identification of drum boiler dynamics. I. P. Eykhoff (Ed.), Identification and System Parameter Estimation, part I, North-Holland, Amsterdam. pp. 78-108. Faurre, P. and J.P. Marmorat (1969). Un algorithme de realisation stochastique. Comptes Rendus Acad. Sci. (Paris) Ser. A, 268, pp. 978-981. Faurre,~L. (1970). Identification par minimisation d'une representation Markovienne de processus aleatoire. Lecture Notes on Mathematics, 132, pp. 85-107; Springer Verlag, Berlin:Faurre, P.L. (1976). Stochastic realisation algorithms. In R.K. Mehra and D.G. Lainiotis (Eds.). System Identification: Advances and Case Studies. Academic Press, New York. pp. 1-25.

482

G. A. van Zee and O. H. Bosgra G.A. van Zee and O.H. Bosgra

Gevers, M.R. and T. Kailath (1973). An innovations approach to least-squares estimation - Part VI: Discrete-time innovations representations and recursive estimation. IEEE Trans. Autom. Control, 18, pp. 588-600. Gold~B. and C.M. Rader (1969). Digital processing of signals, McGraw-Hill. New York. Gupta, N.K. and R.K. Mehra (1974). Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations. IEEE Trans. Autom. Control, 19, pp. 772-783. Huber, P.L. (1972). Robust statistics: a review. Ann. Math. Statist., 43, pp. 10411067. -Jaakso, U. and U. Nurges (1975). Comments on "Approximate linear realizations of given dimension via Ho's algorithm". IEEE Trans. Autom. Control, 20, 302. Jenkins, G.M. and D.~ Watts (1968). Spectral analysis and its applications. Holden-Day, San Francisco, Ca. Kailath, T. and R. Geesey (1971). An innovations approach to least-squares estimation, Pt. IV: Recursive estimation given lumped covariance functions. IEEE Trans. Autom. Control, 16, pp. 720-727. Kalman, R.E. (1965). Linear-Stochastic filtering theory - reappraisal and outlook. Proc. Brooklyn Symp. on System Theory, pp. 197-205 . --Kwakernaak, H. and R. Sivan (1972). Linear optimal control systems. WileyInterscience. New York. ChaV' 1.2. Ljung, L. (1978). Convergence analysis of parametric identification methods. IEEE Trans. Autom. Control, 23, 770-783.---Ljung, L. and A.J.M. van Overbeek (1978). Validation of approximate models obtained from prediction error identification. Proc. 7th IFAC World Congress, Helsinki, pp. 1899-1906. Molinari, B.P. (1975). The stabilizing solution of the discrete algebraic Riccati equation. IEEE Trans. Autom. Control, 20, 396-399. Morf, M. (1974). Fast algorithms for multivariable systems. Ph.D. diss., Dept. Electr. Eng., Stanford Univ. Powell, M.J.D. (1964). An efficient method for finding the minimum of a function of several variables without calculating derivatives. Computer Journal, 7, pp. 155162. Rissanen, J. (1971). Recursive identification of linear systems. SIAM J. Control, ~, pp. 420-430. Rissanen, J. and T. Kailath (1972). Partial realization of random systems. Automatica, 8, pp. 389-396. Rissanen, J. (1973). Algorithms for triangular decomposition of block Hankel and Toeplitz matrices with applications to factoring positive matrix polynomials. Math. Comput., 27, pp. 147-154. Rosenbrock, H.H. (1970). State-space and multivariable theory. Nelson. London.

Soderstrom, T. (1974). Convergence properties of the generalised least squares identification method. Automatica, 10, pp. 617-626. -Soderstrom, T. (1975). On the uniqueness of maximum likelihood identification. Automatica, 11, pp. 193-197. Stewart, G.W. (1973). Introduction to matrix computations. Academic Press, New York. Tse, E. and H.L. Weinert (1975). Structure determination and parameter identification for multivariable stochastic linear systems. IEEE Trans. Autom. Control, 20, pp. 603-613. Van Zee, G.A. and O.~ Bosgra (1979). Approximate linear realizations via singular value decomposition. Report 155-N,Delft Univ. Technol., Lab. Meas. Control, Delft, The Netherlands; submitted for publication. Willems, J.C. (1971). Least-squares stationary optimal control and the algebraic Riccati equation. IEEE Trans. Autom. Control, 16, pp. 621-634. Zeiger, H.P. and A.J:-McEwen (1974). Approximate linear realizations of given dimension via Ho's algorithm. IEEE Trans. Autorn. Control, .!.2., pp. 153. APPENDIX I Define the discrete fourier transformation and its inverse as in Gold and Rader (1969), N-I -j2nfk/N { } E z(k) e fE O,I, .• ,N-I , (I) k=O IN-I j2nkf/N } z(k) N E Z(f) e kdO,I, .. ,N-1 ,(2) f=O with z€RZ and Z€CZ . The sample spectrum is defined as in e.g. Jenkins and Watts (1968), Z(£) =

S(f) = ~ Z(f) z*T(f) ,

(3)

which can via (I) be expressed as N-I S(f) E Rsum(k) e-j2nfk/N k=O with

(4)

Rsum(k) = R(k) + R(N-k)

k={1,2, .• ,N-I}

(5)

RSum(O) = R(O) Now S(f) as defined in (3) is positivehermitian, Le.

(6)

*T~

~f S(f)~f~O and real V~IC

Z, Vf €{O,I, .. ,N-I } (7)

Firstly it will be proved that Rsum(k) is positive-real on [O,N-I], Le. N-I ~ ,*, TR~sum(l' _t) 1\' ~ O an d rea 1 rr\ E 1\ r l\ l'€ C , 1 t i,t=O i=O, I, .. ,N-I . (8) Using the inverse transformation of (4), the positive-real condition (8) can be rewritten as N-I N-I N-I E E \~T ~ E S(f) e j2n(i-t)f/N\ t->0 , (9) i=O t=O 1 N f=O

483

The use of realization theory Realization theory in robust identification

N-I N-I N-I APPENDIX I I E E E A~T ej2TIfi/N S(f)A e-j2TIft/N~0, t (10) N f=O i=O t=O 1 If an approximate model (A,B,C,~) is obtained N-I (11 ) for the sample covariance function and if it is .!.. E not positive-real, then a nearest positive-real N f=O Using the positive-hermitian property of S(f) model (A,B,C,R6) can be computed via a in (7) it follows that the inequality (11) is constrained optimization procedure with parameter Ro, objective function J(Ro )' true and consequently Rsum(k) in (8) is positive-real. J(R) = tr{(R*-R ) (R*-R )T}, 00000 Further, if zeros are contained in the and with a constraint function C(Ro ) that keeps observation sequence, such that Rc within the feasible region, where (A,B,C,Rc) z(k)=O, k=N-M, ... ,N-I, M5N/2 is positive-real. If the algorithm is started then with a feasible initial the constraint R(k)=O, k=N-M, ... ,N-I, function can be chosen as a barrier function and by (5) .!..

RJ,

isum(k) = R(k)

k=O,I, •• ,M.

C(Ro )

Noting that Rsumik) is positive-real on [O,M] it follows that R(k) is positive-real on [0 ,M].

FI+u I

I

I c +x I

Fig.

Fig. 2

F 2+u 2

I

=

~ !1-Ai(M)!

+ O(Ro ) ,

preventing M (4.3) from having eigenvalues !Ai(M)!=1 and preserving the positive sign in (4.2). This is done by defining o(R o ) as o(Rc)=O if (4.2) is true o(Rc)=OO if (4.2) is not true. A feasible initial R~ is obtained by increasing the diagonal elements in R~. This optimization problem is well-conditioned because good initial values are available, the parameter set is small, and the computation of J(Ro) and C(Ro) is efficient.

3

stirred tank

a.

covariance function

b.

sample covariance function

Fig. 3

singular values of HE containing the sample correlation function

484

G. A. van Zee and O. H. Bosgra

__

_ _ _ ...

~=-C.

Fig. 5

Fig. 4

a. second order realization b. third order realization

,

h

Y2 Y2

a.

exact

b.

second order realization

c.

adjusted second order realization

d.

third order realization

:j~:1 ---

Y1 U l

.. .

g ,;

enlargements of R

g

0

,;

I

,;

\:- L I

gj o

g

0

lG

,,'0

6~D~

~~

,;

Fig. 6

impulse response a. exact b. second order model c. third order model

9.

10- 6

0 1.

Q.

0 10- 6

0

0 4.

0 10- 7

0 2.49 10 - 5

6 6.96 10Q

0

0

9.66 10- 7

0

0

7 l.95 10 -

].6j 10 - 7 ]. l O 10 - 4

2 6 8.98 10-

Q

0

0

0

0

9 . 86 10- 7

0

0

7 l.94 10 -

l

- ]. 90 10- 1 l.ll 10- 3

Table I

intensity Q of the innovations process

Qe exact Q2 second order model Q3 third order model