A SPECTRAL REPRESENTATION METHOD FOR CONTINUOUS-TIME STOCHASTIC SYSTEM ESTIMATION BASED ON ANALOG DATA RECORDS P. 1'urprasetiol and S.D. Fassois2 Report UM-MEAM-91-16 01991 by P. Nurprasetio and S.D. Fassois All rights reserved. 1 Graduate Research Assistant, Mechanical Engineering and Applied Mechanics 2 Assistant Professor, Mechanical Engineering and Applied Mechanics

-,I etE "-'. I*s~. 41

1 INTRODUCTION After two decades of almost complete dominance of discrete-time system identification and parameter estimation approaches in the engineering theory and practice, the relevance and importance of continuous-time methods based on analog data have been increasingly recognized (Young, 1981; Young and Jakeman, 1980; Unbehauen and Rao, 1987; Wahlberg, 1989; Neumann et al., 1989; Sagara and Zhao, 1990; Unbehauen and Rao, 1990). A survey of the available literature, however, quickly reveals that the overwhelming majority of currently available continuous-time methods is restricted to the deterministic case, and is thus incapable of handling the more general and practically important class of stochastic systemsl, with which additional difficulties are associated (Young, 1981; Wahlberg, 1989). One of the main difficulties in estimating continuous-time stochastic systems from analog data records is due to our inability to accurately compute time derivatives (of various orders) of the observed random signals (Young, 1981; Wahlberg, 1989; Neumann et al., 1989). If that were possible, continuous-time versions of many discrete methods could be constructed; see, for instance, the works of Balakrishnan (1973), Bagchi (1975), and Pham-Dinh-Tuan (1977), who developed and analyzed the Maximum Likelihood method for this case and proved its asymptotic optimality, and also Unbehauen and Rao (1987), and, Priestley (1981), who discussed the autoregressive model case. The use of the so-called State Variable Filters (SVFs), that are extensively used in the deterministic case (Young and Jakeman, 1980; Young, 1981; Unbehauen and Rao, 1987) in alleviating this difficulty is neither trivial nor effective, as their proper selection is not obvious [and in fact it is known to optimally require a-priori system knowledge, in view of which adaptive SVFs have been also suggested (Neumann et al., 1989)], their incorporation affects the stochastic part of the model, and, unlike the deterministic case, cannot prevent the introduction of large errors into the obtained random signal derivatives (Neumann et al., 1989). Attempts to overcome these problems through alternative approaches based on integral (instead of differential) operators have been also considered. Van Schuppen (1983) examined recursive such estimation algorithms for continuous-time autoregressive processes and proved their convergence, while Moore (1988) analyzed the convergence of the continuous-time version of the Recursive Extended Least Squares (RELS) estimation algorithm for AutoRegressive Moving Average with eXogeneous excitation (ARMAX) systems. In order to avoid numerical instability problems due to the pure integration of the observed signals, exponentially stable prefilters have to be, however, used, the selection of which is important for convergence and, also, requires some form of a-priori system information. Another related problem in this case is that of initial conditions, the effects of which 1A large number of engineering and physical systems are, indeed, stochastic in nature. Consider, for instance, the ambient vibrations of a building or structure, the vibrations of a machine tool during cutting, or those of traveling ground vehicles and aircrafts.

do not decay away and can cause significant estimation bias errors (Young, 1981; Sagara and Zhao, 1990). It therefore comes as no surprise that, for almost all practical cases, the stochastic part of a given system is either not estimated at all, or, only a discrete-time representation of it is obtained based on sampled data (Young and Jakeman, 1980; Young, 1981; Neumann et al., 1989; Sagara and Zhao, 1990; Lee and Fassois, 1991). The only, and arithmetically very few, currently available exceptions, appear to be based on an approach that postulates the estimation of continuous-time stochastic models through the estimation of a directly, and approximately, discretized representation by using, once again, uniformly sampled data records (Astrom and Killstrom, 1973; Phadke and Wu, 1974; Pandit and Wu, 1975). The method developed by Pandit and Wu (Pandit, 1973; Wu, 1977; Pandit and Wu, 1983) is, in fact, a very systematic and comprehensive such procedure, that has been successfully used for the solution of a large number of production engineering, random vibration, and other types of problems. A number of additional such methods have been also, though independently, developed in the fields of statistics (Robinson, 1980; Jones, 1981; Solo, 1984) and econometrics (Sargan, 1974; Bergstrom, 1976; 1983), where the assumption of sampled data is practically reasonable. The main advantage of these methods is in their ability to utilize the already available machinery for discrete-time system estimation (Priestley, 1981; Pandit and Wu, 1983; Ljung, 1987). Yet, and apart from being limited to sampled data records, a number of problems are known to be associated with them, namely: (a) The estimates of the continuous-time parameters are asymptotically biased; with the bias errors being due to the approximate nature of the employed discretization [Young (1981); also see the discussion in Ben Mrad and Fassois (1991), and, Lee and Fassois (1991)], and also analyzed in the econometrics literature [Wymer (1972); Sargan (1974); Bergstrom (1976)], (b) additional errors are introduced into the continuous-time parameter estimates due to sensitivity problems associated with the highly nonlinear discrete-to-continuous transformation (Fassois et al., 1990); with the related question of sampling also being recognized as nontrivial (Fassois et al., 1990; Unbehauen and Rao, 1990), (c) the loss of relative order, that is the difference between the degrees of the numerator and denominator polynomials in the system transfer functions, (d) the incorporation of a-priori system information into the estimation procedure is very difficult, if not completely impossible, due to the highly complicated and nonlinear nature of the discrete-tocontinuous transformation (Wu, 1977; Wahlberg, 1989), (e) high computational complexity, which is further aggrevated by the aforementioned nonlinear transformation. In this paper a novel Mazimum Likelihood type method for the estimation of continuous-time stochastic systems from analog data, that utilizes the general Autoregressive Moving Average with eXogeneous excitation (ARMAX) stochastic model structure and block-pulse function (BPF) spectral representations, and overcomes many of the aforementioned difficulties and limitations, is

introduced. The use of BPF signal representations has led to a number of developments in deterministic system theory and identification in recent years [see, for instance, Wang (1982), and, Jiang and Schaufelberger (1985)]; the present work, however, appears to be the first that, building upon them, utilizes BPF expansions in solving stochastic estimation problems. Within the context of this paper the validity of BPF spectral representations for stochastic signals is, first, formally justified, and, then, based on them, as well as a set of linear operations motivated by recent developments in deterministic identification (Jiang and Schaufelberger, 1985), the problem of interest is shown to be transformed into that of estimating the parameters of an induced stochastic difference equation driven by the spectral representations of exogeneous (observable) and endogeneous (unobservable) excitations; the latter essentially being the spectral expansion of a Wiener process. The structural and probabilistic properties of this difference equation are studied, and the endogeneous excitation is shown to be amenable to a first-order Integrated Moving Average (IMA) representation driven by a stationary innovations sequence. Based on this, as well as the structural properties of the stochastic difference equation, an induced, special-form, discrete ARMAX system with parameters that are expressed as linear combinations of those of the original stochastic differential equation, is finally obtained. The mapping between this discrete ARMAX system and the original stochastic differential equation is shown to be a bijective transformation, and, the discrete system stationary and invertible. These features lead to a Maximum Likelihood type estimation scheme that is based on the estimation of the discrete ARMAX system from which the parameters of the stochastic differential equation are subsequently obtained through a simple linear transformation. The linearity of this transformation, which results in reduced computational complexity while also circumventing sensitivity problems and allowing for the incorporation of a-priori information into the estimation procedure, along with the use of analog data without requiring estimates of signal derivatives, and the elimination of direct discretization procedures, are among the main features of this method, the performance characteristics of which are finally evaluated via a number of simulation experiments. The remaining part of this paper is organized as follows: The exact problem statement is presented in Section 2, some preliminary considerations and the validity of BPF spectral expansions for random signals are discussed in Section 3, the induced stochastic difference equation is derived in Section 4, and its structural properties analyzed in 5. The stochastic modeling of the endogeneous signal spectral representation is discussed in Section 6, and the estimation method formulated in Section 7. Simulation results are presented in Section 8, and the conclusions are finally summarized in Section 9.

2 PROBLEM STATEMENT Consider a continuous-time stochastic system with exogeneous (observable) excitation {u(t)} and corresponding response {y(t)} described by the ARMAX stochastic differential equation2: Sc: ___Ey~t" o d'y(t) du(t) -- dk (t) (1)' 1- bdti dt3 dtk i — j=0 0k=O in which the differentiation operations are to be interpreted in the mean-square sense (Jazwinski, 1970), ana 1, cno = 1, and {w(t)} represents an endogeneous (unobservable) zero-mean continuous-time innovations signal with autocovariance function: E[w(s)w(t)] = (o,o)2. 6(s - t) (2) where 6(.) stands for the Dirac delta function. By using the mean-square differential operator D, the system Sc may be equivalently rewritten in the notational form: Sc: a~(D)~ y(t) = b~(D)~ u(t) + cO(D) w(t) (3) with a~(D), b~(D), and c~(D) being the Autoregressive (AR), Exogeneous (X), and Moving-Average (MA) polynomials, respectively, which axe of the forms: D) - D + Dn- +. + aD + ao (4a) bo(D) _ Dnab b Dna n-2 +... blD + bo (4b) cO(D ],:_x Dnal + co na_2 c(D) Dna-,_ +... cOD + cO (4c) The ARMAX system Sc is additionally assumed to satisfy the following standard assumptions: Al. The polynomials a~(D), b~(D), and c~(D) are coprime (irreducibility assumption). A2. The polynomials a~(D) and c~(D) are strictly minimum phase (stationarity and invertibility assumptions, respectively). A3. The signals {w(t)}, {u(t)}, and {y(t)} are Gaussian, with the latter two additionally being continuous in probability and having almost every sample path characterized by finite energy within the observation interval [O,T), namely fj'[x(t)]2 dt < oo almost surely (a.s.) (Jazwinski, 1970). 2The superscript o is used to indicate quantities associated with the actual system and distinguish them from those of corresponding candidate models.

As we will see the Gaussianity assumption is not critical for the development of the estimation method, but only the vehicle in casting the problem into a Maximum Likelihood framework. It may be thus relaxed, and the method formally interpreted within the broader Prediction Error estimation context (Ljung, 1987). The estimation problem considered in this work may be then posed as follows: "Given analog excitation {u(t)} and response {y(t)} data records generated by the system Sc, subject to assumptions A1-A3, over a period of time [O,T), estimate an ARMAX model of the form: Mc: a(D) y(t) = b(D). u(t) + c(D). w(t) E[w(s)w(t)] -2 6(s - t) (5) that matches the system Sc as closely as possible." 3 PRELIMINARY CONSIDERATIONS The ARMAX equation (1) provides a strictly formal representation of the underlying stochastic system, as it is well known that the continuous-time white noise signal is not second-order and its mean-square derivatives fail to exist. A more appropriate, and also convenient for our purposes, representation may be thus obtained by integrating the differential equation (1) na times, which, assuming, for simplicity, zero initial conditions, yields: y(t) + anal y(t')dt' +... + ao... y(t') dt' nc-fold = b,-1 j u(t')dt' +... + bo... u(t')dt' + ~i(t) + + co... j;i(t')dt' (6) na-fold (na-l)-fold In this expression the integration operations are to be interpreted in the mean-square sense as well (Jazwinski, 1970), while {wt(t)} represents the Wiener process: Ad t CV(t) -wfo)ds (7) which, due to the stated properties of {w(t)}, is Gaussian, continuous in probability, with fgT[zi(t)]2 dt < oo (a.s.), and also zero-mean and with autocovariance function (Jazwinski, 1970): E[w(s)zi(t)] = (oy)2. min(s,t) (8) For the development of the estimation method the m-th order block-pulse function (BPF) spectral representations (Kwong and Chen, 1981) of the signals {u(t)}, {wZ(t)}, and {y(t)}, are now

used. Within the observation interval [O, T) these signals may be then expressed as follows3: m u(t) t(t) Ukk(t) = UTxI(t) (9a) k=l y(t) ym(t) - Yk /k(t) = YT@ (t) (9b) k=1l m z(t) -' 7m(t)- E Wk 4k(t) = WyT(t) (9c) k=1 where {Ui}, {Yi} and {WVi} represent the sequences of BPF expansion coefficients of {u(t)}, {y(t)}, and {iv(t)}, respectively, U, Y, and W corresponding vector representations, and x1(t) an m-dimensional BPF vector of the form: (t)[,(t),(t) [... ~,(t)] tT (10) with 44(t) representing the k-th block pulse function of the m-th order BPF function set: 1 for (k- 1)_<t< < kT?k(t) = - (k = 1,2,...,m) (11) 0 otherwise and Tim the block pulse function duration ("width"). For a given signal, say {y(t)}, the BPF expansion coefficients can be computed through the expression: m T kL Yk y(t)(t)dt= m = im y(t)dt (k = 1,.,m) (12) T = T (k-1)T Such BPF expansions are frequently used for the representation of deterministic signals, and it is well-known, that, in the limit, as m approaches infinity, the set (11) is complete, and the signal representation {ym(t)} given by (9b) converges to {y(t)} pointwise for any deterministic squareintegrable signal defined in the interval [0, T) (Kwong and Chen, 1981). Within the context of the present work the question of validity of such expansions for stochastic signals is an apparently legitimate and important one, and is therefore addressed in the following theorem: Theorem 1: The BPF spectral expansions (9a)-(9c) of the stochastic signals {u(t)}, {y(t)}, and {(f(t)} are convergent in the following product measure sense (exemplified for the case of {y(t))): lim (P x A)((w,t): y(t)- ym(t)I > ) } = VE > 0 (13) with w denoting an elementary event of the sample space, P the probability measure, and A the Lebesgue measure on [O,T). O 3Bold-face characters indicate vector/matrix quantities.

Proof: The proof is a direct consequence of the completeness of the set (11) as m -+ oc, the continuity in probability and finite energy for almost every sample path, that is f:[y(t)]2 dt < oo (a.s.), properties of the random signals involved (see assumption A3 and the earlier discussion on {wz(t)} ), and Theorem 2 of Bharucha and Kadota (1970), which states that the expansion of a random process {y(t)}, with t E [O,T), in terms of an arbitrary basis in ~2(T), the space of squareintegrable functions, is always convergent in the above product measure sense for signals satisfying the aforementioned conditions, regardless of the orthogonality of the basis used and the boundedness of the time interval T. 4 AN INDUCED MAPPING AND A STOCHASTIC DIFFERENCE EQUATION REPRESENTATION By substituting the expansions (9) into the system expression (6), and using the operational matrix form of the BPF spectral representation for integration (see Appendix), we arrive at the following system representation: na na —1 na-1 yT E a _iFix,(t) = UT b~_o _Fi+ l(t) + WT Z C1_Filri(t) (14) i=O i=o i=O in which Fk represents the m x m operational matrix form for k-fold integration (see Appendix). By canceling the vector @(t) from both sides of this equation and performing the algebra, we obtain (na + 1) equations E0,..., E, of the form: ((im) I+k f- I+k Ek n (T/m)' ao +k n, (Tm)'+' Ujfi+l +k+l-j) E: (i + 1)! na-i Yjfi,l+k+l-j)= > ((i + 2)! n —i b Ujfi+i,l+k+l-j) ~i=O ~ j=1 ( i=O j=( na- (T/m)i l+k i=O + ((i )! Cna-l —i E fi,+k+1-) (k =O,l,...,na) (15) with I representing an arbitrary positive integer and fk,i a quantity defined in the Appendix [Eq.(A.3)]. In a manner analogous to that used by Jiang and Schaufelberger (1985) for deterministic systems, by performing the linear operation: (-1)( 7k ) Enak (16) on the set of equations {Ek}, we obtain the following stochastic difference equation: ZAt *Y+i = BV~ U,+i + + CZ * W+i (17) i i=O i=O i=

which is independent of any initial conditions, and, through appropriate reindexing, may be expressed as: SD: A~(B)~ Yk = B(B) - Uk + C~(B) * Wk (k = 1,2,,m) (18) In this equation {Yk}, {Uk} and f{Wk} represent the BPF expansions of the response, exogeneous excitation, and Wiener process, respectively, and A~(B), B(B), and C0(B), polynomials in the backshift operator B (B = Yk _ Yk-l) and of the respective forms: A(B) - AO + A B+... + A.Bnal +A.Bna (19a) =\ t - 1 ao L~m~ —l I ao Dna1(19n) B (B) - B o +B - B +. +..Bna- + B.Bna (19b) na - a-1 1 CO(B) Cna + Cna-1 B +..+ C B"a+ Bna (19c) The coefficients {A}i0nao may be shown to be related to the coefficients of the original stochastic differential equation {a}'na-l through the expressions: na Ao (T/m)a3 / A? = a j(T +l j l!anj (-1) ( k ) fif/, na-k-i+l) (i 1, 2,... n) (20b) j=o I-k=0 and similar expressions relate the rest of the coefficients in (19) with those of the original equation Sc. By defining the continuous and discrete parameter vectors as: A o oa +F a-[ana al a ]T [(na + 1) x 1] (21a) b-~[b~_l... bo bo]T [na x 1] (21b) c [C na-1 C1 [na x 1] (21c) and: A-[An~..... AA n;]T [(na + 1) x 1] (22a) B= [B... B1 B0]T [(na + 1) x 1] (22b)

[C... ]T [(n ) x 1] (22c) respectively, the relations between the discrete and continuous system parameters may be compactly expressed as: A = DA a B = DB b C = DC c (23) with DA being a square matrix with elements determined from (20), and DB and DC submatrices of DA formed by expressing DA in terms of its column vectors as: DA= [ dA(1) I... I dA(na + 1) ] [(n+) x(na+l)] (24a) and defining: DB [ dA(2) I..I dA(n+l 1) ] [(n+l) x na (24b) D= [ dA(1) I *.. dA(na) ] [(n + 1) x na] (24c) These mapping relationships between the continuous and corresponding discrete parameters are summarized in Table 1 for up to fourth-order systems. Equations (23) define a mapping relationship T between the sets: C a {(a, b c) E Rna+l x a xna I with a, b, c of the form (21) witha~ = ca -= 1} (25a) and: D_ {(A, B, C) E 7Zna+l xRn+lXcRna+l I A = DA a; B = DB b; C = DC C; for all (a, b, a) E C} (25b) The nature of this mapping is of particular importance for our developments and is therefore examined in the following lemma: Lemma 1: The mapping T: C -D is a bijective (one-to-one and onto) transformation. Proof: T is a transformation since for any triple (a,b,c) E C expressions (23) define a unique triple (A, B, C) E D. T is also onto since by construction of V every triple (A, B, C) E D is the image of at least one triple (a, b, c) E C. That T is one-to-one follows from the full rank property of the matrix DA; a fact that may be shown by using expressions developed by Kraus and Schaufelberger (1990) in such a way as to rewrite DA as the product of two matrices that may be verified to be nonsingular. Therefore the (sub)matrices DB and Dc are also full rank, and the one-to-one property follows.

5 STRUCTURAL PROPERTIES OF THE STOCHASTIC DIFFERENCE EQUATION S$ Before the induced stochastic difference equation SD of Eq.(18) can be estimated, its structural and probabilistic properties need to be determined. Specifically, issues such as the stationarity of SD, its identifiability, and the first and second order properties of the sequence I{Wk}, are all essential for constructing a proper estimation algorithm. It is emphasized that the study of the last issue is indeed necessary, as the noise dynamics of SD are not completely determined by the polynomial C~(B) alone, but also depend upon the correlation structure of {Wk}. With these ideas in mind we proceed to examine the structural properties of the stochastic difference equation SD first. The following theorem discusses the phase characteristics of the polynomials A0(B) and C&(B) of (18), and therefore the latter's stationarity and "partial" invertibility properties: Theorem 2: Consider the continuous-time ARMAX system subject to assumptions Al and A2, and its corresponding stochastic difference equation SD given by (18). Each one of the polynomials A~(B) and C~(B) of the latter will then be: 1. Strictly minimum phase, provided that its continuous-time counterpart is strictly minimum phase. 2. Minimum phase, provided that its continuous-time counterpart is minimum phase. 3. Nonminimum phase, provided that its continuous-time counterpart is nonminimum phase. Proof: (a) Let us examine the A~(B) polynomial first. (al) Consider the system y(t) = go(D) * u(t) with the strictly proper [see (4)] transfer function g0(D) = b0(D)/a0(D). For a strictly minimum phase a0(D) the system go(D) is asymptotically stable (Chen, 1984), which implies that for every bounded excitation the response will be also bounded, that is: V {u(t)}: IIu(t)1loo < oo = Ily(t)Ioo < 00 (26) where llu(t)lloo suptEJl Ju(t)l with J1 [0, oo). Now consider the corresponding discrete system G0(B) - Bo (B)/Ao(B) induced by the BPF expansions of fixed width T/m. Since G~(B) is (by Lemma 1) unique, its response {Yk} to any given excitation {Uk} will be also unique. An arbitrary bounded excitation {Uk} can be, however, constructed from a bounded ({u(t)} through (12); select for instance u(t) = Uk for (k - l). < t < k.c. From (26) the response {y(t)} to the

excitation {u(t)} will be bounded, and, therefore, the (uniquely-determined) response (Yk} will be also bounded because: m (00 Ik-1 1Jo IlY(t)lloo'llk(t)lloo dt - - IIy(t)l10o J(k-l) dt = - IY(t)loo < ~ (27) where IIYJlloo suptEJ2 IYkl with J2 -[1, 00)4. As a consequence, an arbitrary bounded excitation {Uk} results in a bounded response {Yk}, the system Yk = G0(B) * Uk hence is asymptotically stable, and the polynomial A~(B) strictly minimum phase (Chen, 1984). (a2) For a minimum phase a"(D) the system g0(D) bo(D)/a0(D) is stable in the sense of Lyapunov, and, therefore, its impulse response function {g(t)} bounded: Ilg(t)l1oo < oo (28) The corresponding discrete system G0(B) = B0(B)/A0(B) induced by the BPF expansions has impulse response function: Gk = T g(t) * k(t) *dt (k = 1, 2,..) (29) By taking the norms of both sides of (29), and using (28), we have: IIGklloo < T! 19(t)lloo - [!Ok(t)loo * dt k T T Il(t)lloo (k- dt = 9(t)lloo < 00 (30) which implies that {Gk} is also bounded, and, thus, GO(B) stable in the sense of Lyapunov. As a consequence A~(B) is minimum phase (Chen, 1984). (a3) Now assume that a~(D) is nonminimum phase. Then the system g~(D) is not stable in the sense of Lyapunov, and its impulse response function grows unbounded (1lg(t)llI = cc). By using the continuity of g(t) implied by the fact that g~(D) is strictly proper [see (4a), (4b)], the unboundedness of {g(t)} implies that: V M > 0 3 at least one interval A such that A [(k - 1) k-) for some value of k, and forwhich Ig(t)l > M V t E A (31) 4Notice that although IIXklIo, and Ilx(t)lloo designate different types of norms, the former may be interpreted within the context of the latter by defining x'(t) = Xk for (k - 1) __ t < kT, so that: IIXkIloo = II/(t)lloo.

The impulse response {Gk} of the discrete system G~(B) induced by the BPF expansions will then be, for that particular value of k: Gk =- - g(t) j dt = M (t) > t E A (32) G k k= j () T dt < T M * dt = -M (g(t) < O t E ) (33) Based on this we conclude that: VM>O 3 at least one k such that: IGk > M (34) which implies that {Gk} grows unbounded. As a consequence G0(B) is not stable in the sense of Lyapunov, and A~(B) is nonminimum phase. (b) The fact that the above results hold for the polynomial C0(B) as well may be shown as follows: Rewrite the ARMAX system expression SC [Eq.(6)] in terms of the Wiener process {tiC(t)} in the following notational form: Se~ a~(D)~ y(t) = b~(D)- u(t) + c~(D). t(t) (35) with: O(D) D c~(D) = Cn-lDna + c O2Dna-l +. +cD2 + cOD + 0 (36) By defining: c- Cna-1'''C ~)lO] (37) one may readily show that the vectors c and C are related through the transformation expression: C = DA. (38) which is of exactly the same form as the first of (23) that relates A and a. As a consequence all previous results pertaining to A~(B) are fully applicable to the polynomial C~(B) as well. Corollary 1: Lemma 1 and Theorem 2 imply that the converse of the latter holds for all autoregressive and moving average polynomials with coefficients in D. Corollary 2: An immediate consequence of Theorem 2 and assumption A2 pertaining to the strictly minimum phase nature of a~(D) is that the stochastic difference equation SD is stationary. The structure of the polynomial C~(B) is of particular importance in our developments, and therefore further discussed in the following theorem: 13

Theorem 3: The polynomial C0(B) of the stochastic difference equation (18) has a distinct root at B = 1 and can be factored as: CO(B) = (1-B). C~(B) (39) with C~(B) being minimum phase. o ProofThe fact that OC(B) has a root at B = 1 may be shown by using a known (Jiang and Schaufelberger, 1985) property stating that the sum of the At's is proportional to aO, specifically: Z:2O A = (T/m)naao. Because of (38), this is also applicable to C~(B), and therefore: na ECo = C0(B)B=1 (T/m)nacOO (40) i=o since co = 0 [see (37)]. The fact that the root B = 1 is distinct and C0(B) minimum phase, is a consequence of Theorem 2 which implies that C0(B) has to be minimum phase since 0~(D) is such [based on the definition (36) and assumption A2]. [2 6 STOCHASTIC MODELING OF THE DISCRETE ENDOGENEOUS EXCITATION SIGNAL {Wk} For the development of an estimation method for the stochastic difference equation SD [Eq.(18)], the probabilistic properties of the endogeneous excitation signal {Wk}, defined as:,,mrT Wk= T y i(t)'k(t)dt k E [1, m] (41) need to be also analyzed, and a proper stochastic representation of it developed. Due to the linearity of (41) and the Gaussianity assumption A3, {Wk} will be Gaussian, and thus completely characterized by its first and second-order moments. The first-order moment of {Wk} can be immediately verified to be zero: E[Wk] = T E[tv-(t)] * -k(t) ~ dt = 0 Vk E [1,m] (42) since the Wiener process is itself zero-mean. For the calculation of the second-order moment of { Wk} we proceed as follows: E[lWkl M2] = [T ()d ()(t)dtfT - (" rjk-l)v f |-rmin(s, t) ds. dt 14

Consider the cases: (a) Case k = 1: E[IWk] = m2(0)2 t d + ( fw | | s * ds* dt T2 k-1)-Z -1)-Z J2k-1)-Z -1) T I + II (44) where: m2Q )2 J)km [2 m )2 T2 m3-2 T(2 [6k3i - m (k - 1)22k -(k-1 3 (k-1 (k1 ) La - m (45) = [1 ] T(a)2 (45) By symmetry: [2 3 m( w) (46) and thus: E[Wk2] =k-_) -*T (a~o)2 (47) (b) Case k > l: 2joro 2 kT iT E[WkWI] = T2 I(k-1) T(- tdt ds = (I_ 1 T o)2 (48) (c) Case k <1: T2T [ ] T2 j(k-1)I (-1) T -'T k-l) T )- 1) T )'d' ds ~~~~~~~~~mikl(c) Case ) < I' (,)2= (k~ -2).T. a1 r[WkWl] = { 1 sj~ d s ~ dt T2 dk l) T 1-1)T z - m - T ~ (

Clearly, since this autocovariance is not a function of the relative lag k - 1, the sequence {WZk} is nonstationary. In order to further investigate its structure we define a new sequence {Zk} as follows: Zk Wk - Wk.1 (51) {Zk} is, obviously, zero-mean as well, and with autocovariance E[ZkZI] that may be computed as follows: (a) Case k = 1: E[Z'] = E[Wk] - 2E[WkWkl] + E[Wk-,] - (k- 2)(or)2 -2(k-1- 1)I(ao)2+(k1- I ) (2 7 )2 - To o)2 (52) 3 m W M 3m (0 3m (b) Case k- I = 1: E[ZkZk-l] = E[WkWk-l]- E[Wk.1] - E[fWkWk.-2] + E[Wk-lWk-21] (k - 1 - T (a2)2- T ()2- (k -2- 2)-( )2+ (k - 2 - iT - - (aw()2 (53) (c) Case k -1 > 1: E[ZkZI] = E[WkW] - E[Wk-1Wl] - E[WkWI-l] + E[Wk-VlW-l] (1 T)( )2 - _)i(aTo)2_( -1 - T(o)2+(l- 1T 2 = 0 (54) Observe that E[ZkZl] is not a function of k or 1, but depends upon their difference k - 1, and, by using the symmetry property of the autocovariance, we have: 2 2T(o)2 k-l=~ E[ZkZl] = 1-(ao)2 k- 11 = 1 (55) o0 Ik- 11 > 2 This implies that {Zk} is a stationary MA(1) process (Pandit and Wu, 1983), and may be thus modeled as: Zk =- k + OaNk-1 (56) with {Nk} representing a discrete zero-mean Gaussian innovations (uncorrelated) sequence with variance (orv)2. By comparing the autocovariance of the sequence {Zk} to that of a generic MA(1) process given as: f (+)2[1 +1 (2] Ic-I = O r,(k -1) = (a,.V)29o IC - l = 1 (57) 0 I k-I~>2 16

we obtain the following parameters of an invertible MA(1) representation of {Zk}: 01 -0.267949 (58a) (rO)2 -0.622008 (o~,)2 (58b) m These results lead to the following important lemma: Lemma 2: The BPF expansion series {Wk} of the Wiener process {?1(t)} can be modeled as a nonstationary Integrated Moving Average IMA(1,1) process (Priestley, 1981) of the form: (1 - B) Wk = (1 + 0 B). Nk (59) with 00 given by (58a) and {Nk} being a Gaussian zero-mean innovations sequence with variance given by Eq.(58b). 7 THE ESTIMATION METHOD By substituting the IMA(1,1) representation (59) of the discrete endogeneous excitation signal, and also the form of C0(B) given by (39), into the difference equation SD [Eq.(18)], we obtain the following stochastic difference equation: SoD' A(B) ~ Yk = B(B) ~ Uk + CO(B) (1 + OB) Nk Nk,', i.i.d.V(0,(aN)2) (60) with i.i.d. standing for independently identically distributed. Since 0 is a-priori known, we may reexpress SD in terms of the filtered sequences: UkF (1 + 90B)-' *Uk yF (1 + 90B)-1 *Yk (61) and by additionally normalizing the polynomials A~(B), B~(B), and C~(B) by dividing by A'na) we obtain the following normalized stochastic difference equation: A~(B) = YfF = BO~(B). UkF + CO (B). Nk N'k' i.i.d.A/(0, (caVN)2) (62) with: AO'(B) - 1 + AO'B +... + A~o Bna (63a) A [1 A1... A~a]T A/A7a (63b) 17

B~'(B) -Bo + B B +... + B Bna (63c) B'[Bo' B/'... B~' T _ B/A~ (63d) C (B) 1 + B + + CB +. + C_lBn (63e) _CC ]TA c/C0 A - 1C'-...1 C1 Cna1]T -- C/Cna_ - C/Cna (63f) and: A __ AC0 N A Cna- Nk na- Nk (64) Ana na being a zero-mean innovations sequence with variance: (ON,) =(Cn -1 )2 /AA)2 (N)2 (65) The normalized stochastic difference equation Sil [Eq.(62)] with the stationary zero-mean and uncorrelated endogeneous excitation {Nk } can be now identified as a discrete-time ARMAX(na, na, na - 1) model, which is: 1. Normalized by construction (the leading coefficients of the AR and MA polynomials are equal to unity), 2. Stationary, since A~'(B) is strictly minimum phase by virtue of corollary 2, which guarantees that A"(B) is strictly minimum phase, 3. Characterized by a minimum phase MA polynomial CO'(B) by virtue of assumption A2 and Theorem 3, which guarantee that C~(B) is minimum phase, and is therefore identifiable. The proposed estimation method is based on both the identifiability of S~o and the bijective transformation nature of the mapping between the set of all systems of the form Sp'i and that of continuous-time ARMAX systems of the form Sc [see Eq.(1)]. This latter property is formally given by the following theorem: Theorem 4: The mapping Ta between the sets: Ca - {(a, b, c, (u)2) E RnZ+l xlZna xZfna xZ+ I with a, b, c of the form (21) with a~a = C_1 =1 (66a) 18

and: )' - {(A,B,C', (ro)2) E Rna+l x Rn.a+l x na x X + I with A', B', C' generated through (23), (36), and (63), and (ac,)2 by (58b), (65) by all (a, b, c, (a)) E Ca (66b) is a bijective transformation. 0 Proof: First Ta is a transformation as for any 4-tuple (a, b, c, (a~ )2) E Ca expressions (23), (36), (58), (61), (63), and (65) define a unique 4-tuple (A',B', C', (o;,)2) ) D'. Ta is also onto since, by construction of the set VD, every 4-tuple (A',B', C', (aO,)2) E D' is the image of at least one 4-tuple (a,b,c,(ao,)2) E Ca. The fact that Ta is one-to-one may be shown as follows: Assume a 4-tuple (A', B', C', (ai,)2) E Da. A corresponding 4-tuple (a, b, c, (a~)2) E Ca may be then computed by the following sequence of operations: a DA-1 A' (67a) a = a/ag, (67b) A = DA*a (68a) b = AB1-B'r~ A~ (68b) C~'(B) = (1- B). C'(B) (69a) c -AC-1.,. (69b) c = C/C,0_1 (69c) C = DC (70a) T x 0.622008 x ( ) x (o,)2 (70b) 19

where B',. and C' represent arbitrary na-dimensional subvectors of B' and C', respectively, AB, AC, appropriate n, X n, submatrices of DB and DC, respectively, and a [an...a], c - [c... c]T intermediate (unnormalized) parameter vectors. Due to the nature of these expressions and the full rank property of DA, DB, and DC (and thus of AB and Ac), the 4-tuple (a, b, c, (a, )2) thus determined is indeed unique. O The proposed estimation method is thus composed of two main stages: In the first stage a discrete-time ARMAX(na, n,, n - 1) model of the form [compare with Sli given by Eq.(62)]: MsV: A'(B). yF = B'(B). UF + C(B) - E1 (71) with {E} representing the model's one-step-ahead prediction error sequence, is estimated based on the available sequences {UF }1=l and {YF} k=l, whereas, in the second, the parameter estimates of the continuous-time ARMAX system SC are obtained through expressions similar to (67)-(70). Due to inevitable estimation errors, however, the actually obtained model M.v may not belong to the set D' (MD'f D d), in which case it will have no image within the set C' of continuous-time ARMAX systems. This problem may be dealt with by assigning to MVo that continuous-time ARMAX model Mc E C/ whose image in D' is, in some appropriate sense, closest to that of the estimated Mv',. Within the context of this work we have chosen to achieve that by using the Moore-Penrose pseudo-inverse, and the parameters of the continuous-time ARMAX process may be then estimated as follows [compare with (67)-(70)]: For the estimation of the AR parameter vector a~: a-= DA A, (72a) a = ana (72b) where A' represents the estimate of A~' [the coefficients of the polynomial AO' (B)], a [an.. lo]T an intermediate parameter vector, and a = [,an,... alao]T the vector of the final AR parameter estimates [the estimates of the coefficients of the polynomial a~(D)]. 0 For the estimation of the X parameter vector b~: A-=DA a (73a) b = (DTDB) DgB Af a (73b) 20

where A is the vector of the estimated coefficients of A0(B), B' the vector of the estimated coefficients of B0 (B), and b the vector of the final X parameter estimates [the estimates of the coefficients of the polynomial b~(D)]. o For the estimation of the MA parameter vector c~: C'(B) (1 - B) C'(B) (74a) A (DTDc)-_ DTC' (74b) c = DC (74c) where C'(B) is the estimate of C(B)/Cna, C' the vector of the coefficients of C'(B), C'(B) the estimate of C'(B), c an intermediate estimate, and c = [6nf.. cc]T the vector of the final MA parameter estimates [the estimates of the coefficients of the polynomial c~(D)]. For the estimation of spectral height (ac)2: C = Dcc (75a) ^2 m 1 (An 2 U= T 0.622008 Cn, (}N' ) (75b) Summary of the Estimation Method: Based on the above, the proposed estimation method may be summarized as follows: Step 1: Obtain the BPF spectral representations {Uk} =l and {Yk}=l1 of the exogeneous excitation {u(t)} and corresponding response {y(t)} signals, respectively, by using the operation (12). Step 2: Obtain the filtered representations {UF} and {YF} by using expressions (61). Step 3: Fit discrete ARMAX(na, n, n,a - 1) models of the form (71) to the above filtered representations for successively higher values of n, by using Maximum Likelihood estimation. In each case compute the Bayesian Information Criterion (Akaike, 1977): ne log m BIC(ne) = log(aN,)2 + 2 (76) where ns denotes the total number of estimated parameters, (cfN')2 the estimated variance of the discrete innovations {N2}, and m the length of the BPF spectral representations used. The model 21

that yields the smallest BIC is selected as best. Step 4: Obtain the estimates of the continuous-time ARMAX system parameters through expressions (72)-(75). Remarks: The following remarks are finally in order: (a) Existence and Uniqueness of the Estimates a, b, c, and ^?. The existence and uniqueness of the continuous-time ARMAX system parameter estimates is a consequence of Theorem 4 and the full rank property of DB and DC, which, in turn, guarantees (Barnett, 1990) the existence and uniqueness of their Moore-Penrose pseudo-inverses used in Eqs.(73b) and (74b), respectively. (b) Consistency of the Parameter Estimates. As is well-known, the Maximum Likelihood, or, in fact, the Prediction Error, estimator of the parameters of a discrete ARMAX system is, under mild assumptions, consistent (Ljung, 1987), and, therefore, such is the estimator of SDp [Eq.(62)]. Hence, asymptotically (m -+ oo), the estimated model MyD, E D' (in probability), and since the parameters of the continuous-time ARMAX system are obtained as rational functions of the former [see Eqs.(72)-(75)], the application of Slutsky's lemma (Cramer, 1946) implies that the estimator: 9 A [& T T T 2 ]T (77) will be also consistent, that is: o 9o0o as m — + o (78) where 90 represents the true parameter vector and m the length of the BPF spectral representations used. c 8 SIMULATION RESULTS AND DISCUSSION In this section the performance characteristics of the proposed method are evaluated via numerical simulations. For a given continuous-time ARMAX system of the form (3) and specified exogeneous excitation {u(t)}, the response signal {y(t)} is calculated by integrating each one of the following differential equations: a~(D) ~ yl(t) = b~(D)~ u(t) (79a) a~(D) ~ y2(t) = cO(D) ~ w(t) (79b) 22

and superimposing their solutions: y(t) = yl(t) + y2(t) (79c) For the faithful computation of continuous-time responses the integration step At is selected such that 20 < 7/It < 200, with r representing the smallest system period or time constant. Additional care, is, however, needed for the simulation of the stochastic differential equation (79b). Indeed, the continuous-time innovations {w(t)} is approximated by discrete white noise, and the differential equation is integrated by using a fourth-order Runge-Kutta method. In order for the second-order characteristics of the resulting response to match, at the beginning and end of each integration step, those of the theoretical solution, the variance of the discrete white noise is selected equal to: =2 = 3.6- (.a)2/at (80) with (c,)2 representing the desired spectral height of the continuous-time white noise (Riggs and Phillips, 1987). In all cases examined, the exogeneous excitation signals are composed of trains of pulses with amplitudes forming a sequence of Gaussian independently identically distributed (i.i.d.) random variables with zero mean and unit variance, and the BPF expansions are computed from (12) by using Simpson's composite rule (Burden and Faires, 1985). The estimation of a discrete ARMAX model of the form (71) is then based on the computed BPF spectral records, and the estimated model is validated by examining its predictive ability and the characteristics of its residual (the onestep-ahead prediction error) sequence computed from Eq.(71) for the estimated parameter values. For a good model the latter must be uncorrelated, and this is judged by examining whether its normalized sample autocorrelation lies within the 95% confidence interval of ~1.96/y'; (Pandit and Wu, 1983). Once an estimated discrete-time model is succesfully validated and accepted as an accurate system representation, the continuous-time ARMAX parameters are obtained through the expressions (72)-(75). Estimation accuracy is finally judged in terms of parametric error indices of the form: IIe - 0o II EP = o x 100% (81) with 0 representing a selected parameter vector and 1 i Euclidean norm. Estimation Results The estimation of the underdamped ARMAX(2,1,1) system (System A): (D2 + 2D + 16) y(t) = (D + 10). X(t) + (D + 9) w(t) (82)

is considered first, based on data records that are 100 secs long (T = 100 sees) and generated with white noise sequences having spectral heights (or',)2 = 0.005 and (a')2 = 0.01. The integration step and BPF width were selected equal to T/m = 10At and At = 0.01 sees, respectively. In both of the considered cases discrete ARMAX(2,2,1) models were estimated as statistically adequate, and, as the results of Figure 1, that depicts the normalized sample autocorrelation of the residuals lying within the 95% confidence interval of ~1.96//'iiT, indicate, were successfully validated. From these models the continuous-time ARMAX system parameters were obtained, and the final estimation results are summarized in Table 2, with Figure 2 also comparing the estimated frequency response characteristics of both the b(D)/a(D) and c(D)/a(D) transfer functions to their theoretical counterparts. As these results demonstrate excellent accuracy is achieved, with parametric percentage errors confined to reasonably small values. Similar remarks may be made from Table 3, that presents the results of a Monte Carlo analysis of the method based on twenty data records generated by different seed numbers and (o~)2 = 0.005. Next, the estimation of the overdamped ARMAX(2,1,1) system (System B) described by the equation: (D2 + 3D + 2) *y(t) = z(t) + (D + 10) w(t) (83) is considered based on data records that are 100 secs long and generated with white noise sequences having spectral heights (ao )2 0.00005 and (ai)2 = 0.0001. The integration step and BPF width were selected as in the previous example, that is T/m = lOAt and At = 0.01 sees. In both of the considered cases discrete ARMAX(2,2,1) models were estimated and successfully validated (Figure 3), from which the final estimation results presented in Table 4 were obtained. The achievable accuracy is very good, and the estimated frequency response characteristics nicely match the corresponding theoretical curves (Figure 4). A Monte Carlo analysis based on twenty data records and (ao)2 = 0.00005 (Table 5) further confirms the excellent characteristics of the method. In this final case the ARMAX(3,2,2) system (System C): (D3 + 11D2 + 424D + 1200) * y(t) = (2D2 + 60D + 800) * x(t) + (D2 + 12D + 225) w(t) (84) is considered based on data records that are 37.5 secs long and generated with white noise sequences having spectral heights (or )2 - 0.005 and (a )2 = 0.01. The integration step and BPF width were selected as T/m = 8At and At = 3.125 x 10-3 secs, respectively. Discrete ARMAX(3,3,2) models were estimated as adequate (Table 6) and successfully validated (Figure 5). The final estimation results, corresponding frequency response characteristics, and Monte Carlo analysis of the method [(o )2 = 0.005], are presented in Table 7, Figure 6, and Table 8, respectively, from which very good accuracy is, once again, observed.

9 CONCLUSION In this paper a novel and effective Maximum Likelihood type method for the estimation of continuoustime stochastic ARMAX systems from analog data records was introduced. This method is based on block-pulse function spectral representations, through which the problem is transformed into that of estimating the parameters of an induced stochastic difference equation subject to endogeneous and exogeneous excitations. The study of the structural and probabilistic properties of this equation was shown to further reduce the problem into that of estimating a special-form and identifiable discrete ARMAX system from spectral data. The method was then based on a number of key properties that this discrete ARMAX system was shown to possess, including stationarity, invertibility, and the bijective transformation nature of its mapping relationship with the original continuous-time system. Among the unique features and advantages of the proposed estimation method, that make it especially attractive for applications, are: (a) The fact that neither estimates of signal derivatives, nor direct discretizations that lead to asymptotic bias errors, are used, (b) no prefilters or a-priori information regarding the system dynamics is required, (c) the data are allowed to be in analog form; a fact that, apart from its obvious significance, also implies that non-uniformly sampled and/or missing data can be also accomodated if necessary, and, (d) the relationship between the discrete and the original continuous-time system parameters is linear, so that sensitivity problems associated with highly nonlinear transformations are circumvented, the computational complexity is alleviated, and a-priori system knowledge can be readily incorporated into the estimation procedure. The effectiveness and excellent performance characteristics of the proposed method were finally verified via a number of numerical simulations. ACKNOWLEDGEMENT The financial support of the Indonesian Second University Development Project and the Whirlpool Corporation is gratefully acknowledged.

REFERENCES Akaike, H., 1977, "On Entropy Maximization Principle", in Applications of Statistics, edited by Krishnaiah, P.R., North-Holland Publishing Company, Amsterdam, pp.27-41. Astrdm, K.J., and Kallstr6m, C.G., 1973, "Application of System Identification Techniques to the Determination of Ship Dynamics", in the Proceedings of the 3rd IFAC/IFORS Symposium on Identification and System Parameter Estimation, Eykhoff, P., ed., North Holland, pp.415-424. Bagchi, A., 1975, "Continuous Time System Identification with Unknown Noise Covariance" Automatica, Vol. 11, pp.533-536. Balakrishnan, A.V., 1973, Stochastic Differential Systems, Lecture Notes in Economics and Mathematical Systems, Vol. 84, Springer-Verlag. Barnett, S., 1990, Matrices: Methods and Applications, Clarendon Press, Oxford. Ben Mrad, R., and Fassois, S.D., 1991, "Recursive Identification of Vibrating Structures from Noise-Corrupted Observations, - Part 1: Identification Approaches", ASME Journal of Vibration and Acoustics, in press. Bergstrom, A.R., (ed.), 1976, Statistical Inference in Continuous Time Economic Models, North-Holland/American Elsevier. Bergstrom, A.R., 1983, "Gaussian Estimation of Structural Parameters in Higher Order Continuous Time Dynamic Models," Econometrica, Vol. 51, pp.117-152. Bharucha, B.H., and Kadota, T.T., 1970, "On the Representation of Continuous Parameter Processes by a Sequence of Random Variables," IEEE Transactions on Information Theory, Vol. IT-16, pp.139-141. Burden, R.L., and Faires, J.D., 1985, Numerical Analysis, 3rd ed., PWS Publishers, Boston. Chen, C.T., 1984, Linear System Theory and Design, Holt, Rinehart and Winston, New York. Cramer, H., 1946, Mathematical Methods of Statistics, Princeton University Press. Fassois, S.D., Eman, K.F., and Wu, S.M., 1990, "Sensitivity Analysis of the Discrete-toContinuous Dynamic System Transformation," ASME Journal of Dynamic Systems, Measurement, and Control, Vol. 112, pp.1-9. Jazwinski, A.H., 1970, Stochastic Processes and Filtering Theory, Academic Press. Jiang, Z.H., and Schaufelberger, W., 1985, "A New Algorithm for Single-Input Single-Output System Identification via Block-Pulse Functions," International Journal of Systems Science, Vol. 16, pp.1559-1571. Jones, R.H., 1981, "Fitting a Continuous Time Autoregression to Discrete Data," in Applied Time Series Analysis II, D.F. Findley, ed., Academic Press, New York, pp.651-682. Kraus, F., and Schaufelberger, W., 1990, "Identification with Block Pulse Functions, Modulating Functions, and Differential Operators," International Journal of Control, Vol. 51, pp.931-942.

Kwong, C.P., and Chen, C.F., 1981, "The Convergence Properties of Block-Pulse Series," International Journal of Systems Science, Vol. 12, pp.745-751. Lee, J.E., and Fassois, S.D., 1991, "Suboptimum Maximum Likelihood of Structural Parameters from Multiple Excitation Vibration Data", ASME Journal of Vibration and Acoustics, in press. Ljung, L., 1987, System Identification - Theory for the User, Prentice Hall, Englewood Cliffs, N.J. Moore, J.B., 1988, "Convergence of Continuous Time Stochastic ELS Parameter Estimation," Stochastic Processes and Their Applications, Vol. 27, pp.195-215. Neumann, D., Isermann, R., and Nold, S., 1989, "Comparisons of Some Parameter Estimation Methods for Continuous-Time Models," Proceedings of the 8th IFAC/IFORS Symposium on Identification and System Parameter Estimation, H.F. Chen, ed., Pergamon Press, pp.851-856. Pandit, S.M., 1973, "Data-Dependent Systems: Modeling, Analysis, and Optimal Control via Time-Series", Ph.D. Dissertation, University of Wisconsin at Madison. Pandit, S.M., and Wu, S.M., 1975, "Unique Estimates of the Parameters of a Continuous Stationary Stochastic Process," Biometrika, Vol. 62, pp.497-501. Pandit, S.M., and Wu, S.M., 1983, Time Series and System Analysis with Applications, John Wiley and Sons, New York. Phadke, M.S., and Wu, S.M., 1974, "Modeling of Continuous Stochastic Processes from Discrete Observations with Applications to Sunspot Data," Journal of the American Statistical Association, Vol. 69, pp.325-329. Pham-Dinh-Tuan, 1977, "Estimation of Parameters of A Continuous Time Gaussian Stationary Process with Rational Spectral Density," Biometrika, Vol. 64., pp.385-399. Priestley, M.B., 1981, Spectral Analysis and Time Series, Academic Press. Riggs, T.L., Jr., and Phillips, C.L., 1987, "'Modeling Continuous Noise Sources in Digital Simulations," Simulation, Vol. 48, pp.11-18. Robinson, P.M., 1980, "Continuous Model Fitting from Discrete Data," in Directions in Time Series, D. Brillinger, G.C. Tiao, eds., Institute of Mathematical Statistics. Sagara, S., and Zhao, Z.Y., 1990, "Numerical Integration Approach to On-line Identification of Continuous-Time Systems," Automatica, Vol. 26, pp.63-74. Sargan, J.D., 1974, "Some Discrete Approximations to Continuous Time Stochastic Models," Journal of the Royal Statistical Society, Series B, Vol. 36, pp.74-90. Solo, V., 1984, "Some Aspects of Continuous-Discrete Time Series Modelling," in Time Series Analysis of Irregularly Observed Data, E. Parzen, ed., Springer Verlag, pp.325-345. Unbehauen, H., and Rao, G.P., 1987, Identification of Continuous Systems, North-Holland, Amsterdam. Unbehauen, H., and Rao, G.P., 1990, "Continuous-time Approaches to System Identification -

A Survey," Automatica, Vol. 26, pp.23-35. Van Schuppen, J.H., 1983, "Convergence Results for Continuous-Time Adaptive Stochastic Filtering Algorithms," Journal of Mathematical Analysis and Applications, Vol. 96, pp.209-225. Wahlberg, B., 1989, "On the Identification of Continous Time Dynamical Systems," Proceedings of the 8th IFAC/IFORS Symposium on Identification and System Parameter Estimation, H.F. Chen, ed., Pergamon Press, pp.435-440. Wang, C.H., 1982, "Generalized Block-Pulse Operational Matrices and Their Applications to Operational Calculus," International Journal of Control, Vol. 36, pp.67-76. Wu, S.M., 1977, "Dynamic Data System: A New Modeling Technique", ASME Journal of Engineering for Industry, Vol. 99, pp.708-714. Wymer, C.R., 1972, "Econometric Estimation of Stochastic Differential Equation Systems," Econometrica, Vol. 40, pp.565-577. Young, P., 1981, "Parameter Estimation for Continuous-Time Models - A Survey," Automatica, Vol. 17, pp.23-39. Young, P., and Jakeman, A., 1980, "Refined Instrumental Variable Methods of Recursive TimeSeries Analysis - Part III: Extensions," International Journal of Control, Vol. 31, pp.741-764.

APPENDIX: The operational matrix for integration of signals in the BPF representation Consider the m-th order BPF expansion of the signal {y(t)}, as given by expression (9b). The k-fold integral of {y(t)} may be then expressed as: /... J y(t') dt' XTFk~ (t) (A.1) k-fold where: fk,l fk,2 fk,3... fk,m 0 fk,l fk,2.. — fk,m-1 0 0 0... fk,l and: fk,l = 1 (Vk) fk,i = ik+ - 2 (i - I)k+l + (i - 2)k+l (i = 2,3,..., m) (A.3) The matrix Fk is called the k-th order operational matrix for integration (Wang, 1982). Finally note that the approximation in (A.1) is due to the truncation error associated with the m-th order BPF signal representation, and, due to Theorem 1, converges to zero as m — + oo.

Table 1: Relationships between the coefficients of Sc [Eq.(1)] and S$ [Eq.(18)]. 1st-order system A1 = a + 1 Tao Ao = -al + Tao B1 = bo Bo = i bo C1 = Co Co = -co 2nd-order system A2 = a2 + ~ +:(L)2ao 2ma 6(m) A1 = -2a2 + 3 ( T )2a -- a l(T)2 Ao = a2o- 1 A... + ao B2 -= 1bl + ()2bo B1 = C ( -,1)2bo Bo = 4b )2bo C1 = -2Cl C0 = c, - _Tco 3rd-order system A3= a3 + a2 + 16) a- +24 0()3ao A2 = -3a3 1 Ta2 + (L)2a, + 4 +)3ao A = 3a3- 2 2+( 2 l- 24- o -2ai 1T b)3 B|2 = 1-3J b2 + () + ( )) bo B2 4h-od 2 +- +,)bl + 24(m) b0 0A - 1 T +a + 1+ ((T)2b. 11+ 6(mbo C2 -,=C1 + 2 2)4O bo BC1 =3C22 _C m ( 2C Bo =- p )2bi + 1 (T)3bo _ = - 1 km o 24 C3 = C2 + 2 mC+ -(')2C (i = 3C2 - ITel!(.Z)2CO A _ = -4... 1 + 2 CO ( 2)3 1 + 4th-order system A4- a4+ -1a3 + (_a2 + 4(;)3 a, + 1-Y( ao A0 =-4a4 - T a3 +! (T)2a2 + 5 T 3, 13 4 Bn) b + #iT b+1 0 Al = -434 - m + - ()2a2 - ( T)3al -C 6m) ao Ao2 = -12 T =:a — 4C_3 + a2 + 2(- 60 a) __=____ _1T- 5 T2 3 13 ()4bo Tb3 12 60 T Aoa - T a3 4- + a -)2b- 1T3b. + _ - 7_5)4bo B4 = T__b3 -2- _ b+T C =C- 3 C + 30

Table 2: Estimation results for System A at two different noise powers. Process Estimated Parameters a Parameters O2 = 0.005 a2 = 0.01 a2 1 1.0000 1.0000 al 2 2.0406 2.0256 ao 16 16.1395 16.0939 bl1 1.0537 1.0623 bo 10 10.1251 10.1007 c1 1 1.0000 1.0000 co 9 8.7383 8.7138 ErA() -o) 0.8992 0.6026 Ev(%) 1.3544 1.1783. EC')- 2.8901 3.1605 0 2 - 0.0112 0.0223 aFor the simulation r./At 157. Table 3: Monte Carlo results for System A. Process Estimated Parameters Parameters Mean Value Standard Deviation al 2 2.0249 0.0833 ao 16 16.0771 0.2961 bl 1 1.0364 0.0293 bo 10 10.0765 0.2234 co 9 8.4713 0.5563 EA(S) -- 0.5025 Ed (%)- 0.8426 EC'(%) - 5.8746 - 2 [ 0.005 0.01144 0.00066 31

Table 4: Estimation results for System B at two different noise powers. Process Estimated Parameters a Parameters i = 0.00005 o: = 0.0001 a2 1 1.0000 1.0000 a1 3 3.0570 3.0618 ao 2 2.0181 2.0194 b, 0 0.0014 0.0021 bo 1 1.0203 1.0217 ci 1 1.0000 1.0000 co 10 9.6879 9.7009 EA ( %) 1.5995 1.7324 El(%) - 2.0381 2.1774 EZC'(%) 3.1054 2.9757 oI 7[2 -- 1.1084 x 10-4 2.2154 x 10-4 aFor the simulation r/At = 50. Table 5: Monte Carlo results for System B. Process Estimated Parameters Parameters Mean Value Standard Deviation al 3 2.94573 0.12459 a0 2 2.00387 0.09451 bl 0 0.00029 0.00312 bo 1 1.00906 0.01792 co 10 9.24614 0.62962 BEr(,) 1.50911... EEr"%) — 0.90632 Cr.(%..) 7.53856 |.a2 [5 x 10-5 | 1.14245 x 10-4 6.92103 x 10-6 || Table 6: Order determination results for System C (a2 = 0.005). System order n' | BIC 2 -8.16820 3 -8.77149 4 -8.77039 32

Table 7: Estimation results for System C at two different noise powers. Process Estimated Parameters a Parameters'2 0.005 2 = 0.01 a3 1 1.0000 1.0000 a2 11 11.7901 11.9525 al 424 433.4398 434.9380 ao 1200 1268.5162 1304.6710 b2 2 2.0402 2.0266 b1 60 64.3514 65.2344 bo 800 822.3662 827.7895 Co 1 1.0000 1.0000 cl 12 12.7478 12.8618 co 225 223.1065 224.7139 EA(%) 5.4345 8.2691 E.(%) - 2.8402 3.5249 E () 0.9035 0.4030 a2 - 0.0103 0.0205 aFor the simulation r/At ~ 100. Table 8: Monte Carlo results for System C. Process Estimated Parameters Parameters Mean Value Standard Deviation a2 11 11.3989 0.2347 al 424 431.7177 4.7444 ao 1200 1225.5114 67.7825 b2 2 2.0780 0.0281 bl 60 62.0788 1.5323 bo 800 824.2760 23.5585 cl 12 12.6146 1.1445 co 225 232.4035 13.7799 EAr(o) 2.0944 Ept (0/o) 3.0371 E (S) - 3.2970 w "o'' 0.005 0.01023 0.00045 33

LIST OF FIGURES Figure 1: The normalized sample autocorrelation of the discrete residual sequence (System A). Figure 2: Frequency response curves of the estimated continuous-time system (System A). Figure 3: The normalized sample autocorrelation of the discrete residual sequence (System B). Figure 4: Frequency response curves of the estimated continuous-time system (System B). Figure 5: The normalized sample autocorrelation of the discrete residual sequence (System C). Figure 6: Frequency response curves of the estimated continuous-time system (System C).

0.4 0.4 0.3 0.3 ~ 0.2 o 0.2 AS 0.1 _c 0.1, -- I 0-~ I -0.2 -0.2 -0.3 -0.3 -0.4 -0.4 - 0 10 20 30 40 0 10 20 30 40 L A G LAG (a) aw 0.005 (b) aTv = 0.01 Figure 1: The normalized sample autocorrelation of the discrete residual sequence (System A).

5.0 20.0 - theoretical - theoretical - - - -estimated, 2 = 0.005 0 - - -estimated, 2 = 0.005 00,X 00atd0.o o0.0 ey s/\timated; OrW = 0.01. - - - estimated, tw = 0.01 -20.0 -5.0 -'o' * * *,, _-40.0 o -i.o.0. 5o - 60.0 -15.0 A-100.0 -20.0 -120.0 -25.0 -140.0 0.0 5.0 10.0 15. 20 0.0. 1.0 15.0 20.0 FREQUENCY [ rod/sec] FREQUENCY [ rod/sec] (a) Transfer Function b(D)/a(D) 5.0 20.0 theoretical - theoretical - - - -estimated, e2 = 0.005 - - - -estimated, ao2 = 0.005 0.0 estimated, w= 0.01 00 - estimated, = 0.01 o0.o~ s}o lo~o ls~o 20.0 o~o s~o. s-20.0.0 -5.0 -40.0o ~ -10.0 " -60.0 ~~~~~~~~~~~~~-20.02~~~~~~~~-150~~~~.0 ~-100o -25.0 -120.0 0.0 5.0 10.0 15.0 20.0 0.0 5.0 10.0 15.0 20.0 FREQUENCY [rod/sec] FREQUENCY (rod/sec] (b) Transfer Function c(D)/a(D) Figure 2: Frequency response curves of the estimated continuous-time system (System A).

0.4 0.4' 0.3 4' 0.3 z Z I _ 0.2 o 0.2.. —I.J A, 0.1 _. 0.1 ~ 0.0 H I I 0.0 -0.1 I I*. I _ -NC -0.2 _-. 2 -0.3 -0.3 -0.4 A 0.4 0 10 20 30 40 0 10 20 30 40 LAG LAG (a) aw = 0.00005 (b) a, = 0.0001 Figure 3: The normalized sample autocorrelation of the discrete residual sequence (System B).

0.0 50.0 - theoretical - theoretical -estimated, av 0.00005 - — estimated, 2g = 0.00005 -20.0 |-\ - estimated, a =..w = 0.0 - - - estimated, a 00001 -20.0 w 0. 0 0.0001 -40.0 r -50.0 cm -60.0 14 -100.0 <.. -80.0 -150.0 -100.0. -200.0 0.0 20.0 40.0 60.0 50.0 100.0 0.0 20.0 40.0 60.0 80.0 10.0o FREQUENCY [rod/scc] FREQUENCY [ rod/sec] (a) Transfer Function b(D)/a(D) 20.0.20.0 theoretical theoretical - - -estimated, 2 - 0.0005 0.0 - - - -estimated, a2 = 0.00005 10.0 T T 10.0 i- - - estimated, w = 0.0001 -20.0 0.0.0' -40.0 -iO. -60.0 -2800 20.0 Q-. -100.0 -30.0 -120.0 -40.0 -140.0 0.0 20.0 40.0 60.0 50.0 100.0 0.0 20.0 40.0 60.0 80.0 100.0 FREQUENCY (rod/sec] FREQUENCY [rod/sec] (b) Transfer Function c(D)/a(D) Figure 4: Frequency response curves of the estimated continuous-time system (System B).

0.4 0.4 i 0.2 _ 0.2 m 0.1 _ 0.1 0.0' 0.0.0 a o0. ~ I a,-0I -.2Figure 5: The normalized sample autocorrelation of the discrete residual sequence.2 (System C). C~ ~ (se C ~.

0.0 20.0 theoretical theoretical -estimated, o2 - 0.005 - - estimated, -2 0.005 -5.0 - - - estimated, 0.01. timat, 0.01 -20.0 -10.0 a- -5..5. —.0 2.-0.0 -20.0 -100.0 FREQUENCY [rod/sec] FREQUENCY [ rod/sec] (a) Transfer Function b(D)/a(D) -- 10.0 20.0 theoretical - theoretical —...estimated, 2- 0.005 —..estimated, u2 0.005.- - - estiated ~=:o o.o 01 -0.0 -etimt-a ted =o,~-30~~~~~~~~~~~.0 -20.0 Ldw, -40.0 a. -60.0 -30.0 -00.0 -100.0 0.0 10.0 20.0 30.0 40.0 50.0 0.0 10.0 20.0 30.0 40.0 50.0 FREQUENCY [rod/sc] FREQUENCY [(rod/sec] (b) Transfer Function c(D)/a(D) Figure 6: Frequency response curves of the estimated continuous-time system (System C).

UNIVERSITY OF MICHIGAN 111111111 l 22 IlII11111111 3 9015 02229 1614