A Suboptimum Maximum Likelihood Approach to Parametric Signal Analysis

S. D. Fassois Assistant Professor, Department of Mechanical Engineering and Applied Mechanics, The University of Michigan, Ann Arbor, Ml 48109-2125

K. F. Eman Associate Professor, Mechanical and Nuclear Engineering Department, Technological Institute, Northwestern University, Evanston, IL 60201

S.M.Wu Reid and Polly Anderson Professor of Manufacturing Technology, Department of Mechanical Engineering and Applied Mechanics, The University of Michigan Ann Arbor, Ml 48109-2125

1

A Suboptimum Maximum Likelihood Approach to Parametric Signal Analysis A computationally efficient approach to stochastic ARMA modeling of wide-sense stationary signals is proposed. The discrete estimator minimizes a modified version of the likelihood function by using exclusively linear techniques, and thus circumventing the high computational complexity of the Maximum Likelihood (ML) method. The proposed approach is thus easy to implement, requires no explicit second order statistical information, and is shown to produce high quality estimates at a very modest computational cost. A recursive version of the algorithm, suitable for on-line implementation, is also developed, and, modeling strategy issues discussed. The effectiveness of the proposed approach is finally established through numerical simulations and comparisons with other suboptimum schemes.

Introduction

Signals monitored on a process plant often contain information that is useful in characterizing the dynamics of the plant, its operating condition, and possible failures or malfunctions. Signal analysis techniques, dealing with the extraction of information from measured signals, have therefore become increasingly popular in recent years and several approaches and applications have appeared in the literature (Pandit and Wu, 1983; Braun, 1984; Braun, 1986). Despite the fact that nonparametric, frequency-domain methods that produce reasonable results for certain types of signals are still dominant in signal analysis practice (Wu et al, 1980; Kay and Marple, 1981; Romberg et al, 1984), parametric methods have gained significant momentum lately as their advantages over their nonparametric counterparts are being recognized (Wu, 1977; Kay and Marple, 1981; Braun, 1986). Nevertheless, autoregressive (AR) models, which are generally inadequate in fully capturing the signal dynamics, are often used, and this leads to an analysis procedure that only partially exploits the advantages of parametric modeling. The reason for this selection is clearly bifold: It has to do with the availability of several, easy to implement, computationally efficient AR modeling methods, as well as the difficulties associated with the estimation of general autoregressive moving-average (ARMA) processes (Kay and Marple, 1981; Pandit and Wu, 1983). Indeed, the most straightforward approach for estimating the parameters of an ARMA process is based on the Maximum Likelihood (ML) principle according Contributed by the Dynamic Systems and Control Division for publication in the JOURNAL OF DYNAMIC SYSTEMS, MEASUREMENT, AND CONTROL. Manuscript

received by the Dynamic System and Control Division September 1987; revised manuscript received April 28, 1988. Associate Editor: J. Beaman.

to which a likelihood function is formulated first, and an optimization procedure is subsequently invoked for its maximization (Box and Jenkins, 1970). Despite the nice properties of the resulting estimator (consistency, efficiency, asymptotic normality (Astrom and Bohlin, 1965)), this technique is rather complicated and computationally expensive, since the likelihood function is highly nonlinear and an iterative optimization scheme required. Moreover, the likelihood function may have several local maxima, and convergence to any one of those leads to completely erroneous results (Astrom and Soderstrom, 1974). As a consequence, alternative, suboptimum techniques, yielding good estimates at a more manageable computational load, are sought. A rather common approach is the one based on the extended Yule-Walker equations. This method requires the estimation of the autocovariance function of the signal at a relatively high number of lags, and its performance varies significantly from one process to another (Kay and Marple, 1981). Another rather simple method has been proposed by Graupe et al. (1975), and is also used by Pandit and Wu (1983) for the estimation of initial parameter values which are subsequently refined by a nonlinear optimization scheme. An alternative procedure has been proposed by Durbin (1961), and is often referred to as the two-stage least squares method (Kashyap and Nashburg, 1974). This approach has been however found to be inefficient, and a refined version, designed to overcome this difficulty at the expense of an increased computational load, has been more recently proposed by Mayne and Firoozan (1982). In this paper a novel, suboptimum maximum likelihood ARMA modeling method, based on exclusively linear techniques, is introduced with the objective of achieving high

Journal of Dynamic Systems, Measurement, and Control Copyright © 1989 by ASME

JUNE 1989, Vol. 111/153

Downloaded 27 Sep 2011 to 141.212.97.74. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

quality estimates at a modest computational cost. The discrete ARMA modeling problem and the proposed approach are presented in Section 2. A recursive version of the method is described in Section 3, and modeling strategy issues are briefly addressed in Section 4. Simulation results and comparisons with existing techniques are finally presented in Section 5.

w[t]

X[t]

0 (*-')

Time S9fies

°(2 ARMA Process

01 («•')- 1

The Discrete ARMA Modeling Algorithm A discrete-time stationary signal [x[t]} is said to be an ARMA (n, m) process if it can be represented by a stochastic difference equation of the form (Pandit and Wu, 1983):

; . z - ' Q°(Z-1) 4 i+ £*;.*-

1 - r > m a x ( n , m)+m (27) where n,m are the orders of the ARMA model to be fitted. Step 2: Obtain the initial MA parameter estimates {8k0) £L l by solving equation (25). Step 3: Obtain (x,_i [t]} by filtering the signal [x[t]} by the MA polynomial 6,_,(^~') (equation (10)). Step 4: Estimate the AR parameter vector 0,- through (16). Step 5: Estimate the MA parameter vector 6/ through (17). A block diagram representation of the algorithm is depicted in Fig. 2, where flow of information is represented by thick lines and signal flow by thin ones. The following remarks are finally in order: (a) Steps 3 through 5 of the algorithm can be iterated until convergence of the estimated parameters is achieved. It has been however found (Fassois, 1986) that the improvement thus obtained is often marginal, and therefore going through Steps 1 through 5 once will be sufficient in many cases. (b) The computational burden of the proposed algorithm is equal to 2 linear LS operations of order p and n, respectively, 1 linear system of equations of order m, 1 recursive filtering of order m, and 1 set of equations requiring l/2'm'(m— 1) multiplications. On the other hand, the ML approach requires the equivalent of 1 LS and 1 filtering operations per iteration of the nonlinear optimization algorithm; the number of the required iterations being impossible to specify in advance. (c) The square matrices of equations (16) and (22) are sample covariance matrices, and therefore symmetric and positive semidefinite. This implies that square-root algorithms (Graupe et al., 1980), such as Cholesky decomposition or Householder transformations, offering better stability and numerical accuracy at a reduced computational load, may be used for the evaluation of the respective estimators. In addition, the computational load required for the solution of equation (25) may be reduced from 0(/w3) multiplication operations, required by standard algorithms, to 0(m2) operations by taking advantage of the Toeplitz structure of the square matrix (see Section 3).

3

The Recursive Modeling Algorithm

A recursive version of the proposed algorithm, useful for cases where a signal model needs to be estimated on-line, as the process evolves with time, is now presented, along with its Fast Kalman-type version that leads to a further reduction of the associated computational load.

V-m+l

=

3.1 The Recursive Algorithm. The recursive suboptimum ML algorithm consists of the following steps:

-/,.

Step 1: Estimation of the inverse function l[t]. Ir-

'~m J

^ "mo _,

(25) Journal of Dynamic Systems, Measurement, and Control

A recursive least squares version of (22) is used. In order to make the estimator capable of tracking time-varying signals, expression (21) is modified to: JUNE1989, Vol.111/155

Downloaded 27 Sep 2011 to 141.212.97.74. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

IM 4 arg min £ y[t,i]• (eAR[i\)2

(28)

(=1

where y[t,i\ is a forgetting profile defined as:

4 = [0O10O2. . . .0Oi]T ( l < ; < m ) is then recursively solved as follows: • For/=1,2 m — \:

t

x

y[t,i\& I I

^

TM=1

(29)

y=/+l

! \[t]} being a forgetting factor sequence introduced in order to enhance the tracking capability of the algorithm. The recursive inverse function estimator then takes the form: i[t]=i[t-i]+K[t].(x[t]-ttR[t]4[t-i}) [pxl] (30a) R[f-l]'f«M [pxl] (306) K[fl = R[t]= •

(36)

W/= 5 / + 1 - £ £ • ? ,

(37a)

«i = - p - ( / + i ) - a r * e /

(376)

7/

(37c)

-Pi+i-gi'n

~0o,+ (w,/X,M; "0(i+l):

(37d) W,/X;

(37e)

1

w

e,.+ (H,A;)*g/

], u^xp] *

[R[f-1]-

(30c)

where the dimensions of the vector/matrix quantities are indicated in the brackets, and the forgetting profile is updated through the recursion: \[t] =

,.X[/-l] + (l-Xo)

(30d)

where the constants X0,X[0] are typically selected to be slightly smaller than one. For the initialization of (30) the selection I[0] = 0, R[0] = c^lp ( a » 0 ) may be made. Step 2: Estimation of the initial MA vector 60[t]. This is accomplished by solving (25), which may be compactly written as:

LmV\'60[f\ = *M

(31)

with obvious symbol definitions. The solution of this equation using standard techniques requires 0(m3) multiplications. Lm[fl is however characterized by a non-Hermitian, Toeplitz structure, and Zohar's algorithm (Zohar, 1974) can be used in order to reduce the computational load to 3m2 multiplications. This is done by dividing both sides of (31) by / r _ m so that: (32)

0o = d„ where:

(32a) dm(i) 4 -/ r _ m + ,.// r _ m (326) is obtained (the time argument has been suppressed for simplicity). Furthermore the notation: ^mV'n// = *r-'m + i-j/*r-n

Po

P-1

Pi

Po

P-1

• •

• •

P-m+l

P - 1

• •

• •

P-m + 2

'm-l Pm-2 Pm-i

Moi=d/

156/Vol.111,JUNE1989

(37/) 7,/X, (37*)

\+i=\-«;7//^/

where the double hat denotes order reversal for the elements of a given vector, and: T a, 4 [ p _ , p _ 2 (l