A LINEAR MULTI-STAGE METHOD FOR MULTIVARIATE ARMAX PROCESS IDENTIFICATION - PART II: EFFECTIVE STRUCTURE/PARAMETER ESTIMATION AND PERFORMANCE EVALUATION S.D. Fassois1 and J.E. Lee2 Report UM-MEAM-90-05 @1990 by S.D. Fassois and J.E. Lee All rights reserved..4 -,,* 1 Assistant Professor, Mechanical Engineering and Applied Mechanics 2 Graduate Research Assistant, Mechanical Engineering and Applied Mechanics

1q23 Z, w 4";; y 4 ""

A LINEAR MULTI-STAGE METHOD FOR MULTIVARIATE ARMAX PROCESS IDENTIFICATION - PART II: EFFECTIVE STRUCTURE/PARAMETER ESTIMATION AND PERFORMANCE EVALUATION S.D. Fassois* and J.E. Lee** Department of Mechanical Engineering and Applied Mechanics The University of Michigan Ann Arbor, Michigan 48109-2125 ABSTRACT In the first part of this paper (Fassois and Lee, 1990) the main and guaranteed-stability forms of a linear multi-stage method for the estimation of parameters in multivariate ARMAX processes represented in the fully-parametrized pseudo-canonical form were introduced. In this second part the case of unknown model structure is considered, and a novel structure and parameter estimation procedure that avoids numerically sensitive structure estimation techniques, while using exclusively linear operations and overcoming the enormous computational requirements that render alternative approaches unrealistic for many practical problems, developed. The proposed estimation procedure is based on the BIC criterion, the fully-parametrized pseudo-canonical form, and a very fast correlation-type form of the parameter estimator that uses the second-order moments of the measured data as a pseudo sufficient statistic. This yields remarkably low computational and memory storage requirements, since the original data need not be stored, and only a portion of the estimation calculations has to be repeated as the model order is varied in search of a statistically adequate representation. The performance characteristics of the proposed linear multi-stage method are finally evaluated via numerical simulations. * Assistant Professor, Mechanical Engineering and Applied Mechanics ** Research Assistant, Mechanical Engineering and Applied Mechanics 1

1. INTRODUCTION In the first part (Fassois -and Lee, 1990) of this paper the main and guaranteed-stability forms of a linear multi-stage method for the estimation of parameters in multivariate ARMAX processes represented in the fully-parametrized pseudo-canonical form were introduced. The main form of the method is based on the replacement of the original nonquadratic Prediction Error (PE) estimation problem by an appropriate sequence of simpler quadratic PE and/or linear subproblems, for which closed-form solutions can be obtained. By achieving low computational complexity, requiring no initial guess parameter values, and overcoming the local extrema/wrong convergence problems associated with nonlinear (such as the PE) techniques, this linear multi-stage method overcomes some of the main drawbacks associated with multivariate ARMAX process identification. In addition, the guaranteed-stability form of the estimator overcomes the often encountered instability problems of alternative techniques, and allows for the estimation of all types of multivariate ARMAX processes, including those characterized by MA roots close to the unit circle. In this second part of the paper the practically important problem in which both the parameters and structure of an appropriate (or statistically adequate) multivariate ARMAX representation are unknown and need to be estimated from input/output observations, is considered. The problem of structure/parameter estimation is clearly far from being trivial, as, in addition to the difficulties encountered within the context of standard parameter estimation techniques (excessive computational complexity, algorithmic instabilities, local extrema/wrong convergence, and the need for initial guess parameter values), one also faces the problem of estimating the process Kronecker indices (in case of a canonical process representation) or the autoregressive (AR), moving average (MA), and exogeneous (X) orders (in case of a pseudo-canonical representation). In early identification methods the structure and parameter estimation problems were tackled separately, with the former preceding the latter. The appropriate ARMAX structure was determined by first estimating the structure of a Matrix Fraction Description (MFD) of the process based on rank tests of covariance matrices or linear dependence tests on the rows of the estimated process Hankel matrix, and then transforming that into the corresponding ARMAX representation (Tse and Weinert, 1975; Guidorzi, 1981). As is however well-known (Gevers and Wertz, 1987), the estimation of the Kronecker indices based on these approaches is not only time-consuming, but also involves numerically sensitive procedures, and its resulting accuracy is of critical importance for the consistent estimation of the process parameters, as inconsistent estimates are obtained if the indices are not correctly determined (Caines and Rissanen, 1974). In addition, problems may exist with the transformation of the MFD to the ARMAX representation, as it is known that some of the structural properties of the process may be thus lost (Gevers and Wertz, 1987). More recently, however, combined structure and parameter estimation methods that are pri2

marily based on the Maximum Likelihood concept and involve the extension of the identification criterion with a penalty term that is a function of model complexity, have been introduced. Based on such ideas the Final Prediction Error (FPE) and Akaike's Information Criterion (AIC) have been introduced by Akaike (1974), but are known to yield inconsistent estimates of the process structure (Shibata, 1976). This problem has been overcome by the introduction of the BIC criterion (Akaike, 1977), which allows for consistent structure and parameter estimation under mild assumptions (Hannan and Kavalieris, 1984). The solution of the combined structure and parameter estimation problem for multivariate ARMAX processes however requires the successive fitting of a (possibly large) number of models to the input/output observations, is characterized by a tremendous computational complexity, and therefore renders itself unrealistic for many practical applications. This problem has been recognized by a numbers of researchers, and in an effort to alleviate it, a combined approach that uses canonical and generic ARMAX forms has been proposed by Hannan and Kavalieris (1984). This approach uses a relatively simple, but statistically inefficient, scheme inspired by Durbin's method (Durbin, 1960) and the BIC criterion for the estimation of the Kronecker indices, and finally the nonlinear Maximum Likelihood method for obtaining statistically efficient parameter estimates within the selected model. In this part of the paper a novel and computationally very efficient procedure for structure and parameter estimation is proposed. The proposed procedure uses exclusively linear techniques, and in addition to overcoming the well-known difficulties associated with parameter estimation in multivariate processes, it circumvents the numerical sensitivity, possible parameter inconsistency, and the tremendous computational complexity problems of alternative approaches. The proposed method is based on a very fast correlation-type form of the linear multi-stage parameter estimator that uses the second-order moments of the observations as a (pseudo) sufficient statistic, the fully-parametrized pseudo-canonical form that reduces the structure estimation problem into the problem of estimating three (AR, MA, and X) process orders, and the BIC criterion. The computational and memory storage requirements of the proposed procedure are remarkably low, as the original observations need not be stored in the memory, exclusively linear techniques are used, and only a relatively small portion of the total estimation calculations needs to be repeated as the model order is varied in search of a statistically adequate representation. The performance characteristics of the proposed linear multi-stage identification method are finally evaluated via numerical simulations, and issues such as the achievable accuracy with different noise-to-signal (N/S) ratios and data lengths, as well as the effectiveness of the structure estimation procedure, examined. This part of the paper is organized as follows: The structure and parameter estimation procedure for multivariate ARMAX processes is developed in Section 2, the performance characteristics of the identification method examined in Section 3, and the conclusions finally summarized in Section 4. 3

2. EFFICIENT STRUCTURE AND PARAMETER ESTIMATION As already mentioned the proposed structure and parameter estimation procedure is based on a very fast correlation-type form of the linear multi-stage parameter estimation method introduced in Part I (Fassois and Lee, 1990) and the BIC criterion in conjunction with the fully-parametrized pseudo-canonical form. The fast correlation-type form of the parameter estimator offers tremendously reduced computational and memory storage requirements, and leads to a computationally very efficient structure and parameter estimation procedure that does not need to store the original data or repeat the full estimation calculations as the model order is varied. 2.1 Exact Problem Statement Recapitulating from the first part of the paper, the structure and parameter estimation problem of concern may be stated as follows: "Given a set of N observations of the m-dimensional input {x[t]} and s-dimensional noisecorrupted output {y[t]} vectors assumed to be generated by the multivariate ARMAX system S:1 S: A~(B) ~ y[t] = B~(B) * x[t] + CO(B) w[t] (1) identify a model of the form: M: A(B, 0). y[t] = B(B, 0) x[t] + C(B, 0) e[t/8] E{e[t/0] eT[t/O]} = lg(0) (2) that best fits the data." In the system of Eq.(l) {w[t]} represents an unobservable s-dimensional zero-mean white noise (innovations) process with covariance matrix Eo, AO(B),B0(B), and C~(B) polynomial matrices in the backshift operator B (defined such that B. y[t] = y[t - 1]), that are refered to as the Autoregressive, Exogeneous, and Moving Average matrices, respectively, and are of the form: AO(B) Is + A(1) ~ B + + AO(na) ~ Bna [s] (3) BO(B) _ Bo(1). B + + B(nb). Bnb [s m] (4) C~(B) Is + C(1) ~ B +. + C(nc) Bn [s x s] (5) where IS stands for the identity matrix of order s, A0(j), B0(j), C0(j) are coefficient matrices, and the quantities in the brackets indicate vector/matrix dimensions. The nonnegative integers na, nb, and nc represent the autoregressive (AR), exogeneous (X), and moving average (MA) orders, respectively. In addition the system S is assumed to satisfy the following standard assumptions: Al. The sequence {w[t]} is an s-dimensional, zero-mean and uncorrelated (white) process that is wide-sense stationary and second-order ergodic with a positive-definite covariance matrix E".'The superscript o is used to indicate the actual system and distinguish it from any given model. 4

A2. The observable input {x[t]} is an m-dimensional, zero-mean, wide-sense stationary and secondorder ergodic process with a positive-definite covariance matrix, and is persistently exciting of sufficiently high order. A3. The input {x[t]} and the white noise {w[t]} are mutually independent. A4. All zeros of the determinant of A0(B) lie outside the unit circle ( "stationarity assumption"). A5. All zeros of the determinant of C0(B) lie outside the unit circle ("invertibility assumption"). A6. The Smith form of [A0(B), B(B),CO(B)] is [I,O 0], or, equivalently, the greatest common left factor of A0(B), B(B) and C0(B) is a unimodular matrix (i.e. a matrix with nonzero constant determinant). A7. Rank[A~(na) B~(nb) C~(nc)] = s A8. The input and output processes {x[t]} and {y[t]} are jointly wide-sense stationary and secondorder ergodic. The last assumption A8 is introduced presently as it is needed in the development of the structure and parameter estimation procedure. The model structure M of Eq.(2) is in the fully-parametrized pseudo-canonical form, which implies that the polynomial matrices A(B, 0), B(B, 0), C(B, 0), and the covariance E(0) are full. In these expressions 0 denotes the 4-tuple = (A, B, C, E) that includes all parameter matrices to be estimated, and e[t/0] the one-step-ahead prediction error corresponding to the model structure. 2.2 The Structure and Parameter Estimation Procedure The structure and parameter estimation procedure consists of the following stages: Stage 1: Calculation of the data covariance matrix. The joint input/output observation vector z[t] E R(s+m)xl is introduced as: z A [yT[t] xT[It] = [yl[tl y[t*...y y'[t] l[t 2[t] * *m[t] ] (6) and its sample autocovariance function R,(k,zj) E R(S+m)X(s+m) computed as follows for 1 < j, k < p: N ZN(k-j) z[t]. zT[t + k - j] Rzz(k - j) (k > j) Rzz(kj) I-=) (7) = N j t=1 z[t + j - k]. T[t](k j) (k < j) This function is used as a (pseudo) sufficient statistic, and the original input/output observations need not be retained for the subsequent stages of the estimation procedure. 0 5

For each particular model (characterized by fixed AR, MA, and X orders) the following Stages 2 through 5 are performed, and the model yielding the smallest value of the BIC criterion (calculated in Stage 5) is selected as statistically adequate: Stage 2: Estimation of a truncated ARX representation. By introducing the matrix 7-2 E Rsxp(S+m) and the vector u2[t] E Rp(s+m)xl as: X/2 1 - Hy() H(1): - Hy(2) Hx(2):.. Hy(p) Hx(p)] (8) u2[t] [yT[t-1] x-1] y [t —2] xT[t —2].:... yT[t -p] xT[t- p] (9) which are permutations of the quantities W11 and ul[t], respectively, used in Part I of the paper, the truncated ARX model structure: P p Ml' E H,(j). y[t - j] = E Hx(j). x[t - j] + e[t/12] (10) j=o j=1 with Hy(O) = I, and C[t/1-(2] representing the associated one-step-ahead prediction error, which serves as an approximate representation of the original ARMAX structure M, may be equivalently expressed as: M': y[t] =(Is ~ uT[t] ) h2 + e[t/'12] [s X 1] (11) with the parameter vector h2 defined as: h2 - col [-Hy(1) H(1) -Hy(2) Hx(2)... - Hy(p) H(p) -_ [ hT(1) hT(2) *.. hT(s) ]T [sp(s + m)xl (12) and the p(s'+ m)-dimensional vector hT(i) representing the i-th row of the matrix'H2. Then, similarly to Stage 1 of the main form of the method, the problem of estimating a model within the structure M1 according to the Conditional Maximum Likelihood (or Prediction Error in the non-Gaussian case) method may be decomposed into a set of s scalar estimation subproblems that can be independently solved. By decomposing the h2(1) vector as: h2(1) - [ h'[r h hT ]T [p( + m) x 1] (13) with the (s + m)-dimensional vector hTT representing the first row of the matrix [-Hy(i) Hx(i)], the first scalar equation of (11) may be expressed as: p y1[t] = ZzT[t - j]' h' + c1[t/'1t2] (14) j=1

and the parameter vector ha obtained by minimizing the quadratic criterion: 1V(12) = N Z(yl[t/1-2] = h)2 N = (Y[(]-j] [] T[t- j] h) (15) t t j=1 j=1 which leads to the following estimator expressions: 1 p 1 N E>Ez[t- k].zT[t- j] ~ h. = Z Ez[t - k] yl[t] (k = 1,.,p) (16) j=l t t,, By letting t vary from -oo to +oo, assuming that z[t]- 0 outside the observation interval (that is for t < 0 and t > N), and using the sample autocovariance function of the joint input/output observation vector z[t] defined in (7), as well as the sample cross-covariance between z[t] and yl[t] defined as: N-k rzyl(k)= N Zz[t]yi[t+k] [(s+ m)xl] (17) t=1 the above estimator expressions may be rewritten as: p ZRzz(k-j) h =h rzy,(k) (k= 1,, p) (18) j=1 or in matrix form: RzZ(O) RZZ(-I) RZZ(-p+l) h rzy, (1) RZZ() RZZ(O) Rzz(-p+2) h/ rzy1 (2) LR,,(p-1) RR(p-2).. RZZ(O) h4 rzyJ (p) Rzz h2(1) = rzyl (19) where the matrix Rzz is block symmetric and Toeplitz since Rzz(-k) = RT(k), and rzy, a vector consisting of the sample cross-covariances at lags k = 1,2,.,p between z[t] and yl[t]. By considering all scalar component equations of the model structure M1, and following procedures similar to the above, the following set of estimator equations may be derived: Rzz h2(i) = rzyi (i = 1,p,s) (20) where rzyi is a vector consisting of the sample cross-covariances at lags k = 1,2, -, p between z[t] and y[t]. By defining as RTy(k) the sample cross-covariance function between z[t] and y[t], and decomposing it as: Rzy(k1) [rz(k) rz2(k) rzy,(k)] [(s + m) x s] (21) 7

the set of estimator expressions given by (20) may be compactly rewritten as: -Hy (1) RZz(0) Rzz(-l)... Rzz(-p + ) Rzy(i) Rzz(1) Rzz()... Rzz(-p+2) - | Rzy(2) HiT(2) (22) Rzz(p - 1) Rz(p- 2)... Rzz(O) R( (p) h RI(P) Based on Eqs.(6), (7), (17) and (21), it may be verified that Rzy(k) consists of the first s columns of Rz(k), so that no calculations beyond those needed for obtaining Rzz(k) for k = 0,1, *,p are required, and the matrix 712 is thus estimated by using the estimator of Eq.(22). 0 Stage 3: Initial estimation of the MA polynomial matrix C~(B). This stage of the procedure is identical to Stage 2 of the main form of the method presented in Section 3 of Part I, or, alternatively, the guaranteed-stability form of Section 4 of Part I, and will not be repeated here. 0 Stage 4: Estimation of the AR and X matrices. For the estimation of the AR and X matrices the model structure M2 defined as 2: nb na M2: YF[t] = EXF[t - j]. col[B(j, 03)] - E YF[t - j] col[A(j, 93)] + E'[t/93] (23) j=1 j=l where 03 represents the vector of all entries of the AR and X coefficient matrices defined as (and being a permutation of 02 used in the main form of the method): 03 - col[ -A(1) B(1)! -A(2) B(2)..-. - A(n) B(n)] (24) c'[t/03] the one-step-ahead prediction error associated with this model structure, and yF[t], XF[t - j], and YF[t - j] the following functions of the input and output vectors: yF[t] - [yT[t] 0 C-1(B) ] col[Is] [s x 1] (25) YF[t - j] _ yT[t- j] 0 C-I(B) [s X s2] (26) XF[t - j] - xT[t - j] ~ C-1(B) [s x m s] (27) with C(B) being the estimate of C~(B) obtained in Stage 3 of the procedure, is rewritten as: M2: YF[t] = U2F[t] 03 + E'[t/03] (28) In this expression the filtered data matrix sequence {U2F[t]} E Rsxsn(s+m) is a permutation of {U1F[t]} of the main form of the method and is defined as: U2F[t] - [YrF[t -1] Xr[t- j1] YF[t- 2] XF[t- 2].. Y[ t- n] XF[t - n]](29) 2For the sake of simplicity in this model structure it has been assumed that n = na = nb. 8

By rewritting the parameter vector 03 as: [v v2_n [sn(s + m) x 1] (30) where vj E ~Rs(s+m)X1 denotes col[ -A(j) B(j) ], and defining the joint input/output filtered observations matrix: ZF[t] - [YF[t] XF[t] [s X s(s + m)] (31) the model structure M2 may be further expressed as: n M2 yF[t] = E ZF[t- j] V + E'[t/03] [s X 1] (32) j=l The parameter vectors vj (j = 1,2, —., n) may be then estimated by using a Prediction Error criterion of the form: J trace( —:E'[t/03]'et[t/03 (33) which leads to the following estimator expressions: NZZZy[t]. ZF[t- k] ZF[t- j]y [t (= 1 ZT[t-, n) (34) =1 t t By letting t vary from -oo to +oo, assuming that the joint input/output filtered observations are zero outside the observation interval, that is ZF[t] O for t < O and t > N, and defining the sample covariance functions of the joint input/output filtered observations matrix ZF[t] as: N (kt= ) F[t]~ ZF[t+k-j] kk - R (k-j) (k > j) Rzz(k,j) = (35) Nt=l Z[t + jk] ZF[t] R (k -j) (k < j) and the sample cross-covariance function between ZF[t] and the filtered output yF[t] as: N-k iz~(k) = ZT[t] yF[t + k] (36) t=1 where R, E S(s+m)xs(s+m), ry E Rs(s+m)X1, and 1 < j, k < n, the above estimator expressions [Eq.(34)] may be put in matrix form as: PRLzz(O) l: zz(-1)... Rzz(-n + 1) v1 izy(i) fRlzz(1) Rzz(O)... Rzz(-n + 2).i2 zy(2) Rzz(n-1) Rzz(n-2) *.. Rzz(O) v n(n) with the leftmost matrix being block symmetric and Toeplitz, since Rzz(-k) = RT(k), as it may be easily verified. The solution of this matrix equation yields the desired AR and X matrix estimates, 9

but before this can be accomplished the sample covariances Rzz(k)}n- and Zy(k)}= need to be computed. Within the context of the proposed procedure this is accomplished as follows: By using the definition of ZF[t] [Eq.(31)], the definitions of YF[t] and XF[t] [Eqs.(26) and (27), respectively], as well as Kronecker product properties (Magnus and Neudecker, 1988), the following expression may be derived: ZF[t] = [YF[t] XF[t]] [ yT[t] C-1(B) xT[t] C-1(B) ] = [yT[t] xT[t]] 0 C-1(B) = zT[t] 0 C1-(B) = [zl[t] z2[t]... Z,+m[t]] C-1(B) (38) in which C(B) represents the estimate of C0(B) obtained in Stage 3 of the algorithm. Based on this, the earlier made assumption of ZF[t] being zero outside the observation interval, and Kronecker product properties, R,,(k) may be expressed as: Rzz(k) = N EZ z[t] ZF[t + k] t zl[t] - ] 0 [C-I(B)]T [zl[t. k]... Zs+m[t + k]] 0 C'-1(B) (39) ZS+m [t] In a similar fashion izy(k) may be expressed as follows: ry(k) = [C-(B)]Tyl[t + k]... ys[t + k]] 0 C- (B) col[I,] t Zs,+m [t] = Rzy(k) col[I,] (40) Based on expressions (39) and (40), the filtered sample covariance sequences {R zz(k)}k=ln+l and {izy(k)}=1 may be computed. Towards this end consider the (u, v)-th block element of Rzz(k), ruv(k) E,X's, which may be expressed as: iuv(k) = N ([1-'(B)]T. Zu[t]).(zv[t + k] C-1(B)) (41) t and by introducing the matrix operator D(B): oo D(B) C-1(B) = ED(i) *Bi (42) i=o with D(O) - Is, Eq.(41) may be written as: riL2(kJ) = z[zDT(i) z~[t - i] Ez[t + k - j] D(j)] t i=O j=O 10

= Z_ Z_ DT (i) D(j) ruv(i + k - j) i=o j=O oo oo oo = ZDT(i). D(i). ruv(k) + E E DT(i) * D(i + j) rv,(k- j) + i=O j=1 i=O oo oo Z Z DT(i + j). D(i). ru(k + j) j=1 i=O D(O). ruv(k) + E b(j) ru(k - j) + E DT(j). ruv(k + j) (43) j=1 j=l where: D(j) = y DT(i) D(i + j) (44) i=O and ru,(k) represents the (u, v)-th scalar element of the sample covariance matrix R,,(k) [Eq.(7)]. Because of the invertibility assumption A5 Eq.(44) may be approximated as: D(j) _ DT(i) D(i + j) (45) i=o with an appropriately selected value of r. By further assuming that: ruv(1); 0 ( 111 >q) (46) for some appropriately selected value of q (an assumption reasonable for a large value of q and stationary processes), the values of the sample covariance {iuv(k)} can be computed. By repeating this procedure for each block element iuv(k) (u, v = 1, *.., s + m ) the whole filtered covariance sequence {R,,(k)} can be computed from the sample covariance R,,(k). Finally, by observing that Rfz(k) in Eq.(40) is included in the first s column blocks of Rtz(k), the required values of i,,(k) are already available and do not need to be separately calculated. 0 Stage 5: Estimation of the MA matrices, the innovations covariance, and the BIC criterion. The MA parameter estimates are updated as in Stage 4 of the basic form of the method discussed in Part I of the paper, and the procedure is therefore not repeated here. On the other hand a procedure for the estimation of the innovations covariance and the BIC criterion based on the second-order moments of the observations is now derived. Towards this end consider the system S of Eq.(1): na nb nc E A~(i) ~ y[t - i] = B~O(i) ~ x[t - i] + E CO(i) ~ w[t - i] (47) i=O i=1 i=O By postmultiplying by yT[t] and taking expectation the following expression is obtained: na nb nc ZA~(i)~ Ry(i) = ZB~(i)~R~y(i) + ZCO(i)~ Rw(i) (48) i=O i=l i=O 11

in which {R,(i)} represents the theoretical autocovariance of the output sequence (R'(i) - E( y[t - i] yT[t] }), {(Ry(i)} the the theoretical cross-covariance between input and output (RoY(i) _ Ef x[t - i] yT[t] }), and {R y(i)} the theoretical cross-covariance between the innovations and the output (Roy (i) E w[t - i] EyT[t] }). Based on the stationarity assumption A4, the output vector y[t] may be also expressed as: y[t] = KP(B)~ x[t] + L~(B)~ w[t] (49) with: 00oo K~(B) = [AO(B)]-1' BS(B) = EK~(j) Bj [s x m] (50) j=l L~(B) - [AO(B)]-l * C~(B) = ~ L~(j) Bi [s x s] (51) j=o where L~(O) = I, and by using (49) and assumption A3, the cross-covariance R'y(i) is written as: Ro,(i) = E E{ w[t - i] wT[t - j] }. LT(j) = E~. L T(i) (52) j=o Based on this, Eq.(48) gives: na nb nc j A~(i)' Ryy(i) = B(i)~ Rxy(i) + E CO(i) o L~T(i) (53) i=o i=1 i=O and by using Kronecker product properties the following relations are obtained: col[A~(i) ~ Ry(i)] = [Is A~(i)] col[Roy(i)] (54) col[B~(i) ~ R~(i)] = [Is 0 B~(i)] col[R~y(i)] (55) col[C~(i) Z~~ L~T(i)] = [L~(i) 0 C~(i)] *col[E~] (56) and (53) rewritten as: na nb nc [Is & A~(i)] * col[RY(i)] = n [Is 0 B~(i)] * col[ROY(i)] + E [L~(i) o C~(i)] col[E0] (57) i=O i=l i=O By replacing the theoretical auto and cross-covariance matrices {ROY(i)}, jao and {R Y(i)}nb by their estimates {R.Y(i)}?ao and {Rxy(i)},b~, respectively, that have been already computed as part of 1{R (i)} given by Eq.(21) (since Ro y (i) = [RjT (i)RET (i)]T), and {L0(i)}o by {L(i)}0 O obtained from the estimator expression: y A(j) L,(i -j) = C(i) i = 0, 1, 2,.*., nc (58) j=o 12

which is derived based on Eq.(51) and in which A(i) O for i > na, and replacing the AR, MA, and X matrices by their respective estimates, the following s2 x s2 dimensional matrix equation is obtained from (57): nc na nb [ [L(i) 0 C(i)]] ~ -l[] [I )] col[Ry(i)]- [Is (i)] col[Rxy(i)] (59) i=O i=O i=l which can be solved for the elements of the innovations covariance. Finally, the BIC criterion is evaluated as: BIC(O) = log detlE(O)l + d x -LN (60) where detl.I denotes determinant, 9 the estimated parameter vector, d the number of estimated process (AR, MA, X) parameters, and N the available number of observations. 0 Remarks: (a). For the evaluation of the estimator expressions (22) and (37) in which the leftmost matrix is block symmetric, the computationally efficient Levinson-Whittle algorithm (Whittle, 1963) may be used. (b). The fast correlation-type parameter estimator of the proposed procedure may be valuable also in applications where the structure is known but the parameters need to be obtained from available observations, and a very fast scheme is required. Indeed, the computational savings (in terms of the required multiplication/division operations) of this over the main form of the estimator are very significant for observation records that are not too short, and for long records (N large), in particular, they approach: (S + m)2 ~ (p+ 1) X 100% (s + m)2.p2 + s(s2 na + s.m. nb)2 + (s + m). s3. nc In addition, of course, the memory storage requirements are also significantly reduced. (c). Finally notice that for long data records (for which the computational complexity problem is most severe) a major portion of the calculations (the only that is proportional to the number of observations) is that associated with Stage 1 of the procedure, which, however, does not need to be repeated as the model order is changed. 13

3. PERFORMANCE EVALUATION In this section the performance of the proposed method is evaluated through a number of numerical experiments. In all cases the process response is generated by using mutually independent, pseudo-random sequences characterized by zero mean and approximately flat spectra as the input and innovations signals, and the estimation is carried out by using the basic forms of the method as needed. The emphasis is on assessing: (a) the achievable estimation accuracy at different noise-tosignal (N/S) ratios, (b) the effects of the selected ARX model order p on accuracy, (c) estimation accuracy with different data lengths, (d) the effects of the iterative parameter estimation procedure on the achievable estimation accuracy, (e) the convergence of the (optional) iterative parameter estimation procedure, and, (f) the effectiveness of the model structure estimation procedure. In this context the signal-to-noise ratio in the i-th output channel [N/S(i)] is defined as the standard deviation of the corrupting noise over that of the uncorrupted response in that channel (see Figure 1): N/S(i) = SD(ni) x 100% (61) SD(vi) and is maintained the same in all output channels. The convergence of the iterative parameter estimation procedure is monitored through the index: I101i - 0(i-)IIII Ji= _ 11 1(- ) (62) where II 1 II denotes octahedric vector norm, 1li the estimate of the theoretical process parameter (AR, MA, and X) vector 98 at iteration i, and convergence is assumed as soon as Ji decreases beyond a prespecified threshold value. Identification accuracy is then judged in terms of the normalized parametric error Ep defined as: Ep 1= 01 - 011i x 100% (63) 1100' i, and also by the trace of the estimated innovations covariance: trace [1] = EEe2[t/o1] (64) i=1 t The finally estimated model is validated by comparing the actual process output with the model-based one-step-ahead predictions within a portion of the data set that was not used in the estimation procedure (cross-validation). In addition, the uncorrelatedness of the estimated innovations (residual error) sequence is examined by using the normalized scalar sample auto and cross-correlations defined as: -p2(k) = t(i[t] - ei)(j[t - k] - j) (65) Pj(k) -/E([t] ) (j[t]-j)(65) 14

with {ei[t]} representing the estimated i-th scalar component of the innovations sequence and Ei its sample mean, and subsequently checking their statistical significance for k $ 0 (Pandit and Wu, 1983). In the numerical experiments that follow, the order p of the truncated ARX model used in the method is selected in the range (2.5 r 5) x na, and the parameter estimator is allowed to iterate until convergence is achieved, with the result corresponding to the minimal value of trace [ ] selected as best. Identification Results Test Case I: The parameter estimation problem for the bivariate ARMAX(2,1,1) process of Table 1 from data records consisting of 1000 samples and corrupted at two (1% and 10%o) N/S ratios is considered. The estimation algorithm was allowed to iterate until convergence was achieved, and as Figure 2 indicates, that happened in the third iteration in the 10% N/S ratio case, with the improvement over the first iteration being minimal. The model thus estimated was subsequently validated, and the normalized auto and crosscorrelations of the model residual errors (shown in Figure 3 for the 10% N/S ratio case) are within the theoretical band that indicates statistically insignificant values. The estimated model-based one-step-ahead predictions for both output signals are compared to the actually observed outputs in Figure 4 for the 10% N/S ratio case, where very good accuracy is observed. The full estimation results for both N/S ratio cases are shown in Table 1, and are in excellent agreement with the theoretical values, with the normalized parametric errors being less than 0.04 % for the autoregressive and exogeneous parameters, and less than of 0.2% for the moving average parameters. The estimated frequency response curves of the eight scalar transfer functions corresponding to this system are compared to their theoretical counterparts in Figure 5, where excellent agreement is once again observed in both the 1% and 10% N/S ratio cases. 0 Test Case II: Parameter estimation of the bivariate ARMAX(2,2,1) process of Table 2 is now considered based on data records of lengths N=250, 500, and 1000, corrupted at 10%o N/S ratio. As the estimation results of the table indicate, the proposed method achieves high accuracy even with the data record consisting of only 250 samples. The maximum normalized parametric error in this case is about 1.3%, and naturally decreases to less than 0.08% as the data length increases to 1000 samples. O Test Case III: In this case the problem of structure and parameter estimation for the bivariate ARMAX(3,2,2) process with theoretical parameters given in Table 3 is considered based on data records consisting of 1000 samples and corrupted at two (l% and 10%a) N/S ratios. The model order was succesfully estimated by successively fitting ARMAX(k,k-l,k-l) models and using the BIC criterion, as Figure 6 indicates for the 10%l N/S ratio case. The full estimation results are 15

shown in Table 3, where very good accuracy is observed, with normalized parametric errors of less than 0.03W% for the autoregressive and exogeneous parameters, and less than 2.3% for the moving average parameters. Since algorithmic instability problems were encountered for k > 3, the guaranteed-stability form (Approximation 1) of the estimator was exclusively used in this test case. 0 Test Case IV: In this final case parameter estimation of the bivariate ARMAX(4,2,3) process of Table 4 is considered based on data records consisting of 1000 samples and corrupted at two different (1% and 10%) N/S ratios. The effects of the iterative parameter estimation procedure on the achievable accuracy in the 10%o N/S ratio case are examined, and the trace of the estimated innovations covariance is shown as a function of the iteration number in Figure 7. From this plot some (small) improvement is observed from the first to the second iteration, but convergence is then achieved. The estimated models are subsequently validated, with the validation results shown in Figures 8 and 9 for the 10% N/S ratio case. Figure 8 shows the normalized auto and cross-correlations of the model residual errors, which deviate from the theoretical band for some lags, but are overall considered to be fine. The model-based one-step-ahead predictions of the two output signals are compared to the actual outputs in Figure 9, where excellent agreement is observed. The complete estimation results for both N/S ratio cases are shown in Table 4, where very small (less than 0.05%) normalized parametric errors are associated with the autoregressive and exogeneous parameters, and errors smaller than 1.2% are associated with the moving average ones. O Remarks: Based on the numerical experiments presented, the following remarks may be made: (a). The proposed method offers high estimation accuracy with both "long" and "short" data records and at different N/S ratios. (b). The truncated ARX model order p used in the method does not seem to significantly affect estimation accuracy as long as a sufficiently high order has been selected. In particular, the rule of thumb used [p - (2.5 - 5) x na] was found to work quite satisfactorily in all test cases examined. (c). The iterative parameter estimation procedure converges fast, and may offer some refinement of the ARMAX model parameter estimates. Although the improvement may be small, this procedure may be a desirable option for cases where maximum accuracy is sought. (d). As expected, instability problems are often encountered during the estimation of multivariate 16

ARMAX models, in particular when higher order models are fitted. The use of the guaranteedstability form of the method therefore seems to be necessary in practice. (e). The structure estimation procedure has been found to be quite effective in determining the correct model order. 4. CONCLUSION In this paper the very important problem of structure and parameter estimation in multivariate ARMAX processes was addressed, and a novel and computationally very efficient estimation procedure proposed. The proposed procedure is based on a very fast correlation-type form of the linear multi-stage parameter estimator that uses the second-order moments of the data as a pseudo sufficient statistic, the BIC criterion, and the fully-parametrized pseudo-canonical form. In addition to overcoming the limitations of the currently available parameter estimation schemes, the proposed procedure avoids numerically sensitive structure estimation techniques, while also offering remarkably low computational and memory storage requirements, as the original data need not be stored and only a portion of the estimation calculations has to be repeated as the model order is varied in search of a statistically adequate representation. The very good performance characteristics of the proposed linear multi-stage identification method were finally verified via numerical simulations with a number of multivariate ARMAX processes, by using short and long data records and different N/S ratios. By offering low computational complexity and guaranteed algorithmic stability, circumventing local extrema/wrong convergence problems, requiring no initial guess parameter values, achieving high accuracy, and providing an effective structure and parameter estimation procedure, the proposed linear multi-stage method overcomes the severe drawbacks that have rendered multivariate ARMAX process identification unrealistic in many practical applications, and provides a pragmatic and effective approach for the solution of this important problem. 17

REFERENCES Akaike, H., 1974, "A New Look at the Statistical Model Identification," IEEE Transactions on Automatic Control, Vol. AC-19, pp. 716-723. Akaike, H., 1977, "On Entropy Maximization Principle," in Applications of Statistics edited by Krishnaiah, P.R., North-Holland Publishing Company, Amsterdam, pp. 27-41. Caines, P.E., and Rissanen, J., 1974, "Maximum Likelihood Estimation of Parameters in Multivariate Gaussian Stochastic Processes", IEEE Transactions on Information Theory, Vol. IF-20, pp. 102-104. Durbin, J., 1960, "Fitting of Time Series Models," Rev. Inst. Stat., Vol. 28, pp. 233-244. Fassois, S. D., and Lee, J. E., 1990, "A Linear Multi-Stage Method for Multivariate ARMAX Process Identification - Part I: The Basic Forms of the Estimator", submitted to the ASME Journal of Dynamic Systems, Measurement, and Control. Gevers, M., and Wertz, V., 1987, "Techniques for the Selection of Identifiable Parametrizations for Multivariable Linear Systems," in Control and Dynamic Systems: Advances in Theory and Applications, edited by Leondes, C.T., Vol. 26, pp. 35-86. Guidorzi, R., 1981, "Invariants and Canonical Forms for Systems Structural and Parametric Identification", Automatica, Vol. 17, pp. 117-133. Hannan, E. J., and Kavalieris, L., 1984, "Multivariable Linear Time Series Models," Advances in Applied Probability, Vol. 16, pp. 492-561. Magnus, J.R., and Neudecker, H., 1988, Matrix Differential Calculus, John Wiley. Pandit, S. M., and Wu, S. M., 1983, Time Series and System Analysis with Applications, JohnWiley and Sons. Shibata, R., 1976, "Selection of the Order of an Autoregressive Model by Akaike's Information Criterion", Biometrika, Vol. 63, pp. 117-126. Tse, E., and Weinert, H.L., 1975, "Structure Determination and Parameter Identification for Multivariable Stochastic Linear Systems", IEEE Transactions on Automatic Control, Vol. AC-20, pp. 603-613. Whittle, P., 1963, "On the Fitting of Multivariate Autoregression and the Approximate Canonical Factorization of a Spectral Density Matrix," Biometrika, Vol. 50, pp. 129-134. 18

PARAMETER PROCESS ESTIMATED PARAMETERS MATRIX PARAMETERS N/S = 1 % N/S = 10 % A(1) 0.4 0.1 0.3999 0.1050 0.4079 0.0927 0.2 0.5 0.1991 0.5046 0.1984 0.4959 A(2) 0.6 0.2 0.6038 0.1942 0.6048 0.1884 0.3 0.4 0.3025 0.3963 0.2966 0.4057 B(1) 1.2 -0.5 1.1968 -0.5071 1.1945 -0.5026 0.6 0.3 0.5976 0.2990 0.6019 0.2936 C(1) 0.6 0.25 0.6104 0.2411 0.6173 0.2291 0.2 0.55 0.1785 0.5629 0.1780 0.5538 ErA(.) _ _- | 0.01051 0.03034 E__(%) 0.00508 0.00198 E,' (%) 0.10685 0.16142 Table 1: Estimation results for a bivariate ARMAX(2,1,1) process at two N/S ratios (Test Case I). PARAMETER PROCESS ESTIMATED PARAMETERS MATRIX PARAMETERS N = 250 N = 500 N = 1000 A(1) 0.4 0.1 0.3885 0.1120 0.3971 0.1013 0.3948 0.1023 0.2 0.5 0.1933 0.4765 0.2050 0.4999 0.2008 0.5019 A(2) 0.6 0.2 0.6047 0.2002 0.5984 0.2016 0.5981 0.1995 0.3 0.4 0.2860 0.3962 0.3051 0.3996 0.2887 0.4034 B(1) -0.244 -0.093 -0.2447 -0.0960 -0.2408 -0.0994 -0.2445 -0.0932 -0.454 -0.470 -0.4541 -0.4686 -0.4505 -0.4657 -0.4510 -0.4706 B(2) 0.502 -0.023 0.5038 -0.0337 0.4990 -0.0265 0.5026 -0.0198 0.416 -0.923 0.4182 -0.9144 0.4169 -0.9325 0.4144 -0.9235 C(1) -0.14 -0.3 -0.1211 -0.1779 -0.1830 -0.3447 -0.1214 -0.3075 0.65 -1.0 0.6384 -1.0434 0.6725 -0.9830 0.6647 -0.9788 EA(%) | 1 0.09950 0.00589 0.00227 E_ _ (%) _ 0.01534 0.01100 0.00011 E_ _ _(%_ ) l l 1.12885 0.42796 0.07119 Table 2: Estimation results for a bivariate ARMAX(2,2,1) process from data records of various lengths (Test Case II; N/S=10%). 19

PARAMETER PROCESS ESTIMATED PARAMETERS MATRIX PARAMETERS N/S = 1 % N/S = 10 % A(1) 0.5 0.0 0.5045 0.0001 0.5100 -0.0014 0.0 -0.5 0.0028 -0.4999 0.0032 -0.4986 A(2) 0.0 0.0 0.0048 0.0014 0.0087 0.0006 -1.0 -0.25 -0.9967 -0.2489 -0.9917 -0.2522 A(3) 0.0 0.0 0.0032 -0.0017 0.0035 -0.0013 -0.5 0.125 -0.4972 0.1236 -0.4954 0.1251 B(1) 0.244 0.093 0.2395 0.0998 0.2394 0.0993 0.454 0.470 0,4537 0.4698 0.4557 0.4693 B(2) -0.502 0.023 -0.4999 0.0230 -0.4987 0.0258 -0.416 0.923 -0.4143 0.9288 -0.4100 0.9228 C(1) 0.4 0.1 0.4145 0.0812 0.4211 0.0804 0.2 0.5 0.1352 0.5001 0.1366 0.5016 C(2) 0.6 0.2 0.5922 0.1333 0.5962 0.1317 0.3 0.4 0.1349 0.3362 0.1415 0.3348 EA(_)___ _0.01075 0.02991 - _ _ ___ 0.00092 0.00477 Es(%) 2.28898 2.18361 Table 3: Estimation results for a bivariate ARMAX(3,2,2) process at two N/S ratios (Test Case III). 20

PARAMETER PROCESS ESTIMATED PARAMETERS MATRIX PARAMETERS N/S = 1 % N/S = 10 % A(1) 0.11 -0.80 0.1114 -0.7996 0.1184 -0.8021 0.03 0.64 0.0304 0.6401 0.0316 0.6371 A(2) -0.25 -1.28 -0.2468 -1.2817 -0.2392 -1.2917 -0.02 -0.11 -0.0177 -0.1109 -0.0116 -0.1164 A(3) 0.02 0.39 0.0208 0.3850 0.0222 0.3686 -0.02 0.01 -0.0192 0.0070 -0.0170 -0.0008 A(4) -0.19 0.05 -0.1905 0.0461 -0.1898 0.0381 0.01 0.04 0.0093 0.0369 0.0082 0.0252 B(1) 0.244 0.093 0.2399 0.0995 0.2400 0.0996 0.454 0.470 0.4541 0.4699 0.4567 0.4692 B(2) -0.502 0.023 -0.5010 0.0239 -0.5007 0.0269 -0.416 0.923 -0.4156 0.9229 -0.4162 0.9211 C(1) 0.2 -0.2 0.2203 -0.2293 0.2254 -0.2324 0.8 -0.6 0.7726 -0.5858 0.7724 -0.5887 C(2) 0.46 0.03 0.4842 -0.0344 0.4910 -0.0359 0.32 0.45 0.3603 0.3634 0.3644 0.3635 C(3) -0.21 -0.16 -0.1939 -0.1716 -0.1902 -0.1739 0.03 -0.32 0.0629 -0.3641 0.0624 -0.3658 EA (_ ) 0.00277 0.04735 Er("%) _ _ 0.00013 0.00161 EC_(%)___ 1.10154 1.19040 Table 4: Estimation results for a bivariate ARMAX(4,2,3) process at two N/S ratios (Test Case IV). 21

LIST OF FIGURES Figure 1: Block diagram representation of a multivariate ARMAX process. Figure 2: Parametric convergence indices as functions of the iteration number (Test Case I, N/S = 10%). Figure 3: Model Validation: Normalized auto and cross-correlations of the scalar residual errors (Test Case I, N/S = 10%). Figure 4: Model Validation: Comparison of the process outputs with their corresponding model-based one-step-ahead predictions (Test Case I, N/S = 10%). Figure 5: Frequency response curves of the bivariate ARMAX(2,1,1) model estimated at 1% and 10% N/S ratios (Test Case I). Figure 6: Estimated BIC values for different model structures (Test Case III, N/S = 10%) Figure 7: Trace[]] as a function of the iteration number (Test Case IV, N/S = 10%). Figure 8: Model Validation: Normalized auto and cross-correlations of the scalar residual errors (Test Case IV, N/S = 10%). Figure 9: Model Validation: Comparison of the process outputs with their corresponding model-based one-step-ahead predictions (Test Case IV, N/S = 10%). 22

w1 2 wt. w INNOVATIONS (COVARIANCE Z ) A (B)C (B) n [t] n[t]..n[t] NOISE x [t] y [t] x [t] 1 y [t] v[t] " r x It] - y I t] EXOGENEOUS UNCORRUPTED CORRUPTED INPUTS OUTPUTS OUTPUTS Figure 1 Block diagram representation of a multivariate ARMAX process.

0.002 -0.001' 0.000 ~ 1 2 3 4 5 6 ITERATION (a) Autoregressive Parameters 0.0010 - 0.0008 0.0006 0.0004 0.0002 0.0000 1 2 3 4 5 6 ITERATION (b) Moving Average Parameters 0.0012 0.0010 0.0008 i';" 0.0006 0.0004 0.0002 0.0000 -T ~ 1 2 3 4 5 6 ITERATION (c) Exogeneous Parameters Figure 2: Parametric convergence indices as functions of the iteration number ( Test Case I, N/S = 10% ).

0.4 0.4 o.2 o.2 o0C 4 - 0.0 0. 0 - 20 0 40- - - -0.4 s I I, -0.4 1 * I.... 0 10 20 30 40 s50 0 10 20 30 40 50 LAG LA G (a) p21 (b) P12 0.4 0.4 0.2 0.2 Z -0.2 -0.0 LAG LAG Figure 3: Model Validation: Normalized auto and cross-co.rrelations. of the scalar:.. Q ~~~~~~~~~~~~o 0.0~~~~esda ros etCs I, NS- 10.0

8.0 4.0 t -: Actual Output Actual Output 8.0 +: Prediction + Prsdiction 2.0 4.0 2.0 o0.0 9o -2.0 -2.0 -4.0 -6.0 -4.0 900 920 040 060 960 1000 go900 920 040 060 980 1000 DATA NUMBER DATA NUMBER (a) output Yl (b) output Y2 Figure 4: Model Validation(: Comparison of the process outputs with their corresponding model-based one-step ahead predictions ( Test Case I, N/S = 10% ).

16.0 50. li O. 10 0 10.0 6C, Q. 50. C: 6.0 6.. I / \ I O~~~~~~~~~~~~CZ I 0 ~ i 0.0 -. -.0 --- C_ —100. Z Z.I C. t.&. - 0.0 () —150. = -200 -6.'6 0 ~1.67 3.14 6.oo 1.67 3.14 FREQUENCY (r/g) FREQUENCY (r/g) (a) Y1/X1 Transfer Function 16.0 60.. 6060. C,jo. iO.~~ -5.0 - -10.0, -— J LAJ 6~~~~~~~~~~~~~10. co -6C.0 _ _ _ _ _ _ -100. 4~~~~~~~~~~~~~~~~~4 -6-.0o,.~-~ 0., -1.0 -) I 50. C, ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~~p -400 1.67 3.14 00 1.7 3.14 FREQUENCY (r/g) FREQUENCY (r/e) (b) Y2/X2 Transfer Function i5. 0 0..... 50., CMI~~~~~~~~~~~~~~~~~~~~CO i0.0 ~- (4-u -400 f 4)~~~~~~ O. 1. 053 —00. -1 6. S. 40'3 L & FREQUENCY Cr/ e) FREQUENCY fr/8) ~~~~~~ ~ ~ ~" — 5 0 R.,,c1 ~. 1o. o — 10.0~~~~~~~~~~~~~~~~~~~~~~~~.' — 1 00 -2'00 6 0I.6731 FREUECY(re)FRQUNC (/e'1~~~~~~~~~~~~~~ -30 43 -cs~g~ oo r~si 5.(4 -sa~d. 00 l~-o FREQUENCY ~~~~~~~~~~r/r) FREQUENCY ~~~~~~~~~~~~~r/~~) — 10.0 ~ ~ d Y/zTrnfe unto

16.0 40.. 20. 10.0 5L O. 6.0 ~'~c3~~~~~~ tw~~~~~L.. -20. C., 0.0 C -40. C., LAJ -6.0 = -io~~~~~~~~g~~~ ~ — 00 — 10.8........*~6~~~~~~~~-0 O.I0 6I.67 3.4-.00 I 61.7 3.14 FREQUENCY (r/r) FREQUENCY Cr/s) (e) Y1/W1 Transfer Function 20.0 100. 0.0.3 ~~~~~~~~ O. Cn w~~~~~~~~~~~ -20.0.-I oo. -4o IIwC= -60.0 o -3o0. I-~~~~~~~~~~~~~~~~~~~~~~L - ~~~~~~~~~~~~~-400 L 0 1.6 oo 37-14006.00 1.57 3.14 FREQUENCY Cr/o) FREQUENCY (r/r) (f) Y1/W2 Transfer Function 100.,... 10.0 6.0 0.. o0 L.J -200.' ~~~~~~~~~~~C.~-:-200. CD - 0.0 - A ~~~~~~~~~~~~' C.,~~~~~~~~~~~~~~~~~~~~~. (/ -300. -0-.0 -400 020 t.00 15 31 FREQUENCYG7 Cr/ 1)FREQUENCY (r/e) FREQUENCY (r/e) (g) Y2/W2 Transfer Function lO.0 20. 0 m 4.0 at ) 1% W an 0%NSranf Fctios etCaeI) N 1% = 00 -20. C.D EC. LaJi — 40. 4C a- I~~~~~~~~~~~~~~~ Z~~~~~~~~~~~~~~~~~~~~~ ILO. 7 3.14 6.00 3714 FREQUENCY (r/o) FREQUENCY (r/o) (h) Y2/W2 Transfer Function at 1% and 10%/ N/S ratios ( Test Case I ). ( ~: actual system, -—' N/S-= 1%, ----- N/S-= 10% )

-2 U -4 (2,1,1) (3,2,2) (43,3) (5,4,4) MODEL ORDER Figure 6 Estimated BIC values for different model structures ( Test Case III, N/S = 10% )

0.0160 0.0158 0.0156 0.0154 0.0152 0.0150 0.0148,. 0 1 2 3 4 5 6 ITERATION Figure 7: Trace [E] as a function of the iteration number ( Test Case IV, N/S = 10% ).

0.4. 0.4 0.2 0.2 (a) Pk (b) P1 ~.4 -: 0.4 r-.. -0.2 - -0.2 - -0.4 -0.4 0 10 20 30 40 50 0 10 20 30 40 50 LA G L A G (c) P21 (d) P22 3C -0.2 -0.2 -0.4 -0.4 0 10 20 30 40 50 0 10 20 30 40 50 LAG LAG (C) P21 (d) P22 rsda e ( T Case IV, N/ = 10% )I

10.0 8.0: Actual Output: Actuol Output +: Prediction + Predlction 4.0 5.0 2.0 - 0.0: I I 0 I.0 9. OE+02 9.2E+02 9.4E+02 9.6E+02 9. E+02 1.OE03 900 920 940 960 980 1000 DATA NUMBER DATA NUMBER (a) output Yl (b) output Y2 Figure 9: Model Validation: Comparison of the process outputs with their corresponding model-based one-step ahead predictions ( Test Case IV, N/S = 10% ).

UNIVERSITY OF MICHIGAN 3 9015 02826 744411111 3 9015 02826 7444