Adaptive estimation using parallel processing techniques

Comput. & Elect. Engng,Vol.2, pp. 203-214.PergamonPress,1975.Printedin GreatBritain.

ADAPTIVE ESTIMATION USING PARALLEL PROCESSING TECHNIQUES* PETER K. S. TAM a n d JOHN B. MOORE Department of Electrical Engineering, University of Newcastle, New South Wales, 2308, Australia (Received 9 August 1974)

Abstract--A parallel processing technique for adaptive estimation is investigated. In particular, for the case of unknown system or statistical parameters denoted by the vector 0 belonging to a finite set {0,, 0_~..... 0N}, the maximum likelihood 0 is determined and denoted 0, and the minimum mean square error state estimate conditioned on this 0, namely ~(tlt, ~) is taken to be the state estimate. Using this approach new estimators are derived which require less computational effort and have less limitations than previous adaptive estimators using parallel processing techniques described in the literature. Results for the case of time-varying unknown parameters are also derived. An example is included of state estimation for a known signal model but with unknown noise statistics. The filter banks are constrained to be time-invariant and so only approximate maximum likelihood parameter estimation is achieved. I. I N T R O D U C T I O N

The Kalman-Bucy filter[l, 2] for the estimation of the states of a linear dynamical system requires an exact knowledge of the system parameters and noise covariances. We consider the adaptive estimation problem of estimating the states when the dynamical and/or statistical model is specified up to a set of unknown parameters, denoted by the vector O. Parallel processing techniques have been applied by a number of authors [3-5] to the adaptive estimation problem and in fact adaptive estimators requiring one hundred or so Kalman filters can be implemented using mini-computers. In essence, the standard Bayesian approach to the adaptive estimation problem is as follows [5]. Assuming that the unknown parameter vector 0 is discrete or suitable quantized to a finite number of grid points {O. . . . . ON},with known or assumed a priori probability for each 0, the conditional mean estimator includes a bank of N Kalman filters where the ith filter is a standard Kalman filter designed on the assumption that 0 = 0i. The filter bank is driven by the noisy signal measurements. The conditional mean state estimate is given by a weighted sum of the states of the Kalman filters. The weighting coefficient of the state of the ith Kalman filter is the a posteriori probability of 0i, which can be updated recursively using the noisy signal measurements and the state of the ith Kalman filter. Unfortunately, for systems with continuous-time measurements, the above Bayesian approach has the drawback that, first, the measurement noise covariance R has to be known, and second, the unknown parameter vector 0 has to be time-invariant. For the special case of systems with discrete measurements, parallel processing estimation techniques have been developed when the measurement noise covariance is unknown, at the expense of either complexity or loss of optimality [4, 6]. In the practically important direction of reducing the complexity of the adaptive estimator, Alspach and Abiri [6, 7] obtain time-invariant state estimators for time-invariant systems with unknown noise covariances by considering a grid of possible time-invariant Kalman gains directly rather than a more involved grid of possible noise covariances in the unknown parameter space. The results can be extended to cases when the noise covariances are time-varying quantities [8]. However, for the estimator of [6] to be close to optimal, it is required that the Kalman filters have reached steady-state and that the number of measurements received be large. Thus, during the transient periods, the estimator operates suboptimally. Some parallel processing estimation techniques are also available for discrete-time systems when 0 is time-varying. Using a combination of digital and analogue techniques, optimal adaptive estimators for the case when the unknown parameter is a scalar Markov sequence of known statistics have been developed [9]. By approximating the a posteriori density of the state vector *Work supported by the Australian Research Grants Committee. 203

204

P.K.S.T.~,M and J. B. MI)I~RI

with a Oaussian probability density, Ackerson and Fu[10] derived a suboptimal estimator for the particular time-varying unknown parameter case in which the input noise or the measurement noise comes from a group of white Gaussian noise sources, which act one at a time, with the transition from one noise source to another being described by a discrete Markov transition matrix. It is clear from the above survey of existing adaptive estimators using parallel processing techniques that it would be worthwhile to investigate any adaptive estimation approach using parallel processing techniques suitable for both continuous-time and discrete-time problems, and which can handle unknown time-varying parameters and unknown measurement noise parameters with but a small addition in complexity over the case of time-invariant unknown parameters and known measurement noise covariance. In this paper, such an approach is investigated. This approach can be briefly described as follows. An unknown parameter vector 0 is defined in such a way that there is one and only one Kalman filter corresponding to a particular grid point 0~ in the 0 space. The states of the Kalman filters are denoted )?(t], 0~). By comparing the relative magnitudes of the likelihood functions of the 0,'s evaluated using the measurement data and £(tlt, 0~) for all i, the most likely 0~ at time t, denoted 0(t), is determined. The state ~;c[t]t, 0(t)] is taken as the estimate of the signal model state x (t). Note that since the likelihood functions of the unknown parameters are used instead of their exact a posteriori probabilities, it is not necessary to assign a priori probabilities to the different O~'s. The advantages of using the approach just outlined can be summarised as follows. First, for continuous-time systems, in contrast to the standard Bayesian approach, our approach does not require that the measurement noise covariance be known exactly. Second, for the discrete-time case, our estimators require less computational effort than alternative known estimators using parallel processing techniques. Third, for both the continuous-time case and the discrete-time case, the various results can be extended in a simple manner to give useful adaptive estimators for the case of time-varying 0. The above advantages are of course useless unless the estimators using the above approach perform satisfactorily. Monte Carlo simulations show that they work very well in minimizing the mean square estimation error when compared to the more complex adaptive estimator for discrete-time problems described in [8]. An outline of the subsequent sections of this paper is as follows. In section I1, we present the discrete-time results for time-invariant unknown parameters. In section III, the results of section II are extended to the case of time-varying unknown parameters. In section IV, the results of sections II, III are extended to systems with continuous-time measurements. II. DISCRETE-TIME RESULTS FOR TIME-INVARIANT UNKNOWN PARAMETERS In this section, we first review appropriate discrete-time optimal adaptive estimator results from[3-5] and then we apply these results to achieve an alternative adaptive estimator which in most applications is simpler and thus more attractive. A simulation example is given to demonstrate the performance characteristics of the alternative estimator. Consider the system

x(tk+,) = cb(tk+t, tkltO )x(tk t+ G(t~ [q,)u(t~ ,

(1)

y(t~+,) = H~tk, ,Ito)x(t~, ,)

t2t

z(tk+,) = y(tk,.)+ v(tk,~)

(3)

where u(.), v(.) are independent zero-mean Gaussian white noise sequences with covariance matrices I and R respectively. The positive definite matrix R will be assumed known in some instances, and not known in others. The initial state x(to) is a Gaussian random vector with mean £o(to) and variance P(tolto, to) and is independent of u(.) and v(.). The entities ck(tk+,, &]to), G(tklto), H(tk+,lto), £(to), P(tolto, to) are completely specified by the parameter vector tO, i.e. they are known if tO is known, possibly unknown if tO is not known. This formulation effectively allows an unknown input noise covariance (although we have assumed E{u(tk)u'(tk)} ~ I) through the intervention of G(tkttO) in (1).

Adaptiveestimationusingparallelprocessingtechniques

205

If 0 and R are specified, the conditional mean state estimate ;(tklt~, 0, R) is of course given from the familiar Kalman filter equations: $(tk*~ltk+,, 0, R) =

~(t~.,lt~, o, R)

=

Y,(t~+llt~, 0, R) + K(tk.,lO, R)f(t~+,]0, R) ~(tk+,ltk, O);(tk Its. 0, R), ~(tolto, 0, R) = ~o(0)

~(t~.,lO, R) = z(tk+~)- H(tk+llO)Yc(t~÷~ltk, 0, R) K(t,,+.lO, R) = P(tk+,lt~, 0, R)H (t~+d0, R)P~ (t~+,[t~, 0, R) t

-1

(4) (5) (6)

(7)

H(tk+,lO)P(t~+,lt~, 41,R)H'(t~+~10) + R

(8)

P(t~.,Itk, 0, R) = &(t~+,, tklO)P(t~ltk, 0, R)&(tk+~, t~10) + t-';(t~10)G'(t~ 10)

(9)

P(tk+,ltk.,, 0, R) = e(t~+l[tk, 0, R) - K(t~+l]~', R)P~(t~+llt~, 0, R)K'(t~+110, R).

(10)

Next we consider that R is known but 0 is unknown. Assuming that ~, the space of admissible values of 0, is discrete or suitably quantized to a finite number of grid points {0, . . . . . ~}, with known or assumed a priori probability p(0~lto, R) for each 0, the conditional mean state estimate £(tkltk, R) is given by [3-5]: P

$(tkltk, R) = ~ 2(t~ltk, 0,, R)p (0, Its, R)

(11)

i=1

P (0, ]tk,R) = c IP~ (tk Ik-,, 0. R)I -''2 exp {-½11~(t~10. R)ll~='.~r ..... .,.-,}P (0, Its-,, R)

(12)

where c is a normalizing constant independent of 0~ so that P

p(O, ltk, R) = 1. i=1

Now consider the more general case where R is unknown with a continuous range R of admissible values. An obvious approach for extending the above results is to first approximate R with a suitable finite set of quantized points {R~. . . . . R,~}. Then with known or assumed a priori probability p (01, RjIto) for each (0, Rs), the conditional mean state estimate would be given from [4] P

M

g(tkltk) = ~ ~ 2(tkltk, 0, Rj)p (0,, R,[tk)

(13)

i=1 j = l

P(O,, Rj[tk) = c'Pz(tkltk-l, 0,. RJ)]-''2 exp {-½ll2(t~ 10. RJ)[I~-'..~ ..... .,.R,,Ip(0,, Rjlt~-d

(14)

where c' is a normalizing constant independent of 01 and Rj. To implement equation (13), P × M Kalman filters are required compared to the requirement of P Kalman filters for the implementation of equation (11). This represents a large increase in the number of parallel processing units when R is unknown. As a first step towards reducing the filter complexity, we replace the above discretization of the Cartesian product • × R with the discretization of an alternative set which we now describe. Through equations (7-10), the product space • × R defines a (continuous) space O of hypotheses of possible Kalman filter configurations, specified by the quadruples {K(.]0, R), q5(.10), .fo(0), H('[0)}. (In many cases, the various time-varying Kalman gains have to be approximated by time-invariant gains to simplify calculations). For each 0 in O, we have the following equations in lieu of equations (4--6). :f (tk+,[t~+., 0) = 2(t~+,lt~, 0) + K(tk+.lo)~(tk+,[O)

;(t~+,ltk, 0) = t~(tk+l, t,[O)2(t~[tk, O)

2(tolto, 0) = go(O)

i(t~+l[O) = z(tk+l) H(t~+,lO):~(tk+llt~, 0). -

(4)' (5)' (6)'

Note that while equations (4--6) can be re-written in terms of 0 as above, it is not possible to

206

P.K.S. TAM and J. B. Moose

do so with equations (7-10). Thus the mapping • x R ~ 6) incurs some loss of information as the specification of 0 alone determines only the Kalman filter configuration but not the error covariance of the state estimate associated with that Kalman filter. It turns out that if O is to be approximated by a set of N grid points {0 . . . . . . 0,,. }, then in many cases N can be much smaller than P × M. As an illustration of the above statement, we point out that in [6], a numerical example for state estimation of a linear time-invariant system with unknown noise covariances is given whereby 2500 grid points in the space of unknown noise covariances can be adequately replaced by 100 grid points in the space of possible steady-state Kalman gains. Using the set of grid points {0, ..... 0~}, the conditional mean state estimate is now given from

£(t~[. ) = ~ 2(tklt~. O, )p(O,]& )

115)

where p(O~]tk) is the a posteriori probability that 0 = 0,. A crude approach for calculating p(O~ltk) is to first observe that the joint a posterion probability that 0 = 0, and R = R~ can be updated recursively using the following equation

P(O,.R~)k) = c"lP:(tklt~ ,. O,.R,)t '/'-exp{-½1l~(t~[o,)ll~.. ' ............... }p(O.R~lt~ ~)

t16)

c" is a normalizing constant independent of O, R~ and P:ltkl& , . 0 ~ , R ) = E{i(tk ]0,)-~'(tk 10,)]Rj } is related to [K (tk ]0, ), Rj] through the following relation (c.f. equations (7), (8)), which holds when p(O, Rj]to)~ 0 where

P~(tkltk ,, 0,. R, ) = [I - H(&]O,)K(tkfO,)] 'R,.

(17)

Then by approximating the Cartesian product space O × R with a set of grid points {(0~, Rj)Ii = 1. . . . . N; j = 1. . . . . M} and applying the theorem of total probability, we can express p(O~lt~) as a marginal probability, viz, M

p(O,l&)= ~ p(O Rjltk). i

(18)

I

This approach of using equations (16-18) however requires considerable on-line calculations as we have to store and update M separate quantities to obtain a single p(O, ltk ). Moreover, a lot of off-l!ne calculations may be required to obtain the a priori probabilities p(O,, Rjlto). An alternative approach is to find an approximate expression for p (0, It~ ). For example, using the techniques described in [6], it can be shown that for scalar measurements and assuming that for each 0~ the range of admissible values of P:(tkltk ~, 0~, R) is the interval (0, ~), then we have the following approximate expression t19~ where c is a normalizing constant independent of 0; and k is a large even number (greater than about 1000). We now propose a simple estimation scheme to yield a suboptimal minimum mean square error estimate for the case when R is unknown with a continuous range of admissible values. Our approach is based on the observation (shown below) that under a simplifying assumption on the range of admissible values of R, then given the measurements up to time tk, it is relatively easy to find [0(tk), 15,(tk)], which we shall denote simply as (0, 15,), that maximizes the likelihood function p(Z~[O, R) of (0, R), where Zk = {z(t.) ..... Z(tk)}, Using the estimate 0 so obtained, our proposed simple estimator is one yielding the estimates 2(<k, 0). (Note that P, is not required other than as a step to find ~J.) Using the fact that E{~(h]O,)~'(tslO,)[R} = [I H(hlO,)K(tj]O,)l 'R, the likelihood function

Adaptiveestimationusingparallelprocessingtechniques

207

p(Z~[O,, R), which is the conditional probability density of Z~ as a function of 0~ and R, can be established as k

l

k

.

2

t

(20) Now we make the simplifying assumption (quite similar to that made in [6] when the measurements are scalars) that for each 0,, R can take any value provided it is a positive definite symmetric matrix. For each 0 = 0, let 1~(0~)denote the value of R which maximizes p (Zk [0~,R). If 1~(0,) exists, it can be obtained by setting [Op(Z~[O~,R)]/OR -~= O, which gives for k > 1,

f{(O,) = ½(ll+ 119

(21)

11= -~l ~ ~(tjlo,)~'(tjlo,)U- I-t(tjlO,)K(tjlo,)]'.

(22)

where k

Substituting the above value of l~(0~) into the expression for p(Zk[O~,R) gives k

p [z, loi, ~(o,)1 = cI-I

-

j=l

H(tjlo,)K(tslO,)[''~}

(23)

where c is a constant independent of [0, 1~(0~)]. Equations (21) and (22), though useful in a number of cases, have to be used with caution because they apply only when 1~(0~) so obtained is positive definite. (Hence they are not applicable for k = 1.) In order to ensure a positive definite l~(0,), we need further assumption on the range of R for each 0 = 0,. Let us assume that for 0 = 0,, R is of the form R = a~R, where R, is a known positive definite matrix and a~ is an unknown scalar which can take any value in (0, oo). Then setting [Op(ZklO,, R)]/0a, = 0 gives

+ ~,Itztt,lo,)ll.,l~(O,) = R_, k [.~_-' -

2

1

(24)

This equation (24) gives a positive definite l~(0~) for all k = 1, 2, 3. . . . Substituting this l~(0i) into the expression for p(Zk]O~,R) again gives equation (23). Note that 1~(0,) and p [Zk]O~,R(O,)] can be updated recursively, and that max {p [Zk ]O,,1~(0,)]} is in fact the global maximum of p (ZklO,, R) with respect to i and R. Our proposed estimator consists of a bank of N Kalman Filters. The conditional innovation processes of the Kalman filters are used to update recursively p[ZklO~, 1~(0,)], by means of equation (23). Then 0 is chosen to maximize p[Z~[O~,1~(0,)], and ~(tkltk, O) is chosen as the state estimate. Such an estimator is sketched diagrammatically in Fig. 1. Under certain circumstances, the computations can be further simplified. For example, if the measurements are scalars, then as maybe seen from equation (23), the problem of maximizing k

p(ZklO, R) with respect to (0, R) reduces to the simpler problem of minimizing { E ~2(tjlO~)[i= j=l

1 ....

N}.

Application to state estimators in unknown noise environment It is interesting to compare the performance of our proposed estimator with the estimator proposed by Alspach and Abiri in [6] for state estimation of a linear time-invariant system when both the plant driving noise covariance Q and the measurement noise covariance R are unknown. In [6], to reduce the complexity of the filter structure, the [Q, R] space is mapped to the space of

208

P. K. S. TAM and J. B. MOORZ

z(k)

If~*er

I£(t

t 8)

!

II

p[4O. R(o)

i ".~ .......

I

I tilter

>p[ZtO #O)

I x(tlt,~) ~ ~ ~

i

•

Select

I

-

---

r~E~r°~,~(@ '

~

] i~'~

* J' / f o r 011 j

s~,~,

....

£(tl t,O)

L

l

f

i

~~tltS)

[ ~lg,ven II h~ i-#~ '

L~f_ ,'N ['~(t J8 ) ' ! k N I

t

I Updot#

[

/

r ~l"lSI°,'~cs~;lr~

Fig. 1, Discrete-tune adaptive estimator for tune-invariant 0. 47

o

~~ .,..,",~AIsoach

E

era/

15

I 2

t 3

t 4

L 5

L 6

1 7

I 8

I 9

I J I I 0 I'112

tt3

1 1 1 I J I i 14 15 16 17 1 8 1 9 2 0 2 ~

Time Fig. 2. Mean square error for state estimation in high noise (G = 0.6, J = 10).

Kalman gains (assumed to be time invariant). The simplifying assumption that the Kalman gains are time invariant leads to sub-optimal performance during the transient period but is still a very useful assumption to make. In order to provide conditions for comparison purposes will also make this assumption for our proposed estimator. Thus we consider the grid {0, . . . . . ON}to be the grid of possible time-invariant (sub-optimal) Kalman gains. We now consider the state estimation of a first-order system where the noise covariances are unknown. The system equations are x(tk+,) = 0.8x(tk) + Gu(tk) z(tk) = x(tk) + Jv(tk )

(25) (26)

where u(.), v(.) are zero-mean independent white Gaussian sequences with unity covariances and G, J are 'unknown' quantities. The unknown Kalman gain for the above system can take a value between 0 and 1. We consider the following simple Kalman gain space consisting of N grid points {01. . . . . ON} as the unknown parameter space. The ith Kalman filter (that is the filter conditioned on 0 = 01) has Kalman gain K(O~)= 0~ where O~(i-0.5)/N.

Adaptive estimation using parallel processing techniques

209

For the simulation results, the signal generating system is first run until steady-state can be assumed. Then measurements z(k) are processed for estimation of x(tk) for k = 0, 1. . . . . 21. Typical performances (the mean-square-error of the state estimates at each tk) are shown in Figs. 2 and 3, using 300 sample paths and 10 grid points in the Kalman gain space. The results show that our estimator actually performs slightly better than that of [6]. To conclude this section we comment that not only does our scheme give better performance than that of [6], at least for the somewhat random examples chosen, but it is certainly considerably less complicated to implement.

~-

0"1

"~ ,'-

/

~

/

Proposed \

~

Ispoch

e;e OZ

estimator

0.05

3 4 5 678

9 I01112131415161718192021 Time

Fig. 3. Mean square error for state estimation in low noise (G = 0.6, J = 0.1).

Ill, D I S C R E T E - T I M E

RESULTS

FOR TIME-VARYING

UNKNOWN

PARAMETERS

In this section, we show that the estimation scheme described in the last section can be extended in a simple but suboptimal manner to handle unknown time-varying (but not slowly varying) parameters and unknown measurement noise parameters. The addition in complexity, which involves feedback from the estimator output to the individual Kalman filters, is very small relative to that required for an optimal minimum square error estimation scheme. For our problem, equations (4'-6') generalize to g[tk+,ltk+,, O(tk*l)]

=

g[tk+,ltk, O(tk)] +

K[tk+,lO(t~+,)]~[t~+~lO(tk)]

g[tk+~[tk,O(tk)] = &[tk+~, tk[O(tk)]2[tk[tk, O(tk)], ~[tk+,lO(t~*,)] =

g[totto,O(to)] = g[O(to)]

z(t~+J- H[tk+~lO(tk+l)].f[tk+,]tk,O(tk)]

(27) (28) (29)

where O(t~) = {0(to) . . . . . O(tk). Here we have to assume that for 0 = 0, the measurement noise covariance R(tk) for v(tk) is of the form R(t~) = a(tk)R,(tk) where a(t~) is an unknown scalar and Ri(tk) is a known matrix. We assume that a(t~), O(tk), O(tk-i) are independent of each other and that a(tk) can take any value in (0, oo). The extimator structure that we are going to propose can be slightly modified to include the case where there is a particular range of admissible values of O(tk-,) for each O(tk). However, simulation results so far indicate that very little improvement in performance can be gained by imposing such restrictions. At t = t~ max {p[ZllO(t, O(to), a(h)]} for each O(tJ can be found in two stages as follows. O(to).a(tD

First, note that we can express p [Z,IO(t,), O(to), a(h)] in an equation similar to that of equation (20) for p[ZklO, R]. Therefore, for each of the grid points {[01(t,), Oj(to)]li, j = l . . . . . N}, p[ZdO(tO, O(to), &(h)]= max p[Z~lO(tO, O(to), a(t,)]} can be found via a pair of equations =(tO

similar to that of equations (23) and (24). Second, O(to) is chosen such that p[Z, lO(h), O(to), &(h)] = max {p[Z~lO(tl), O(h), O(to), a(h)]} where O(to) can in general be different for different O(t)

C.A.E.E., Vol. 2. No. 213--F

210

P.K.S.

'rAM and J. B. MOORE

O(tO. Having thus found 0(to) for each O(tO, the one step ahead prediction state estimate given by the Kalman filter conditioned on each O(t,) is then reset to g[t,[to, 0(to)I, where again .~[tdto, 0(to)] can be different for different O(t~). To recapitulate, based on the 'maximum likelihood' approach as described in the last section and under some simplifying assumptions, we reject, for each O(t,) at t = tt, all the combinations of [0(to), a (tO] except the pair [0(to). c~(L)] which maximizes the likelihood function p [Z, IO(t,). 0(to), ~ (t,)]. In general, at t = t~.. max p[L]O(t,), O(t, ,). ~)(t, ~), c~(t,), ACt, ,)} for each O(tk) can o(tj.

I),~(t~ }

be found in two stages similar to those described above for t - t,. Here (:)(t~ =)={0(to), 0(t0 . . . . . 0(tk =)} and A(t~ ,)={c~(t,) . . . . . a(t~ ~)} and these two sequences can in general be different for different O(h ,). Having thus obtained 0(tk ,) for each O(tk), the one step ahead prediction state estimate given by the Kalman filter conditioned on each O(t~) is then reset to £[tkltk ,, O(t~ ,)]. Resetting the states of the Kalman filters as described above allows us to limit the number of Kalman filters to a finite number N although the number of grid points for O(tk) grows exponentially with k. As an illustration of the above procedure, we shall describe the derivation of the various quantities in more details for the case of scalar measurements (i.e. R(tk)=c~(t~)) and H[tk]O(tk)] = H(tk) independent of O(tk), which gives particularly simple results. At t = t,. we have

p [Z,lO(t,), O(to), a(t,)] = {27ra(t,)} 'r'{1 - H(t,)K[t,lO(t,)]} ''2 exp

f

l[t,lO(t,)]~

I

(30)

where

l[t,[®(t,)] = {1 - H(t,)K [t,l®(t,)]}i2[t d®(t,,)]

(31)

~[t,[O(t,,)] = z(t,) - H(t,)Y,[t,lto, 0(to)].

(32)

Observe that, for each [0(t,), 0(to)l, p = p[ZdO(t,), O(to), a(t,)] is a unimodel function of a(t,). Setting [Op/Oa(t,)] = 0 gives &(t,) = I[tllO(tl)] and p[Z, lO(tO, 0(t,,), 6(t,)] = c]~[t,lO(to)] [. Choosing 0(to) to maximize {lf[tllO(to)]}, the one step ahead prediction state estimate given by the N Kalman filters are then reset to ~f[tllto, 0(to)]. (Note that in this case 0(to) is independent of 0(tl) so that the calculations are much simplified.) By inducticn, we easily obtain for each [0(t~), O(t~ ,)]

d~(tk) = I[tklO(t~), O(tk ,), O(tk 2)]

(33)

p [Z~lO(t~ ), O(t~ ~). (~(tk _~),A (t~)] = c'q~[tk IO(t~ O, O(t~ 2)]1

(34)

where l[tklO(tk), O(tk ,), ~)(tk 2)]:{1--H(tk)K[tk[O(tk)l}y2[tklo(tk ,), ~)(tk ~-)l, 2[tklO(tk ,), O(tk-2)] = z(tk)-H(tk);[tkltk ~, O(tk ,), (9(tk 2)]. The normalising constant c' and the sequence O(tk-2) are independent of [0(tk), O(tk 3]. From equation (34), p[ZklO(tk), O(t~ O, ~)(tk-=), A(tk)] is independent of O(t~). Thus, the problem of finding the global maximum of p [ Z IO(tk), O(tk O, (~(tk 2). a(tk), A(t~ ,)] with respect to [0(tk-,), a(tO] reduces to the problem of finding m!n {22[tklO~(tk ,), e)(tk 2)]}. The one step ahead prediction state estimate given by the N Kalman filters are then reset to ; [tklt~,, (~(tk_,)]. The estimator structure for this problem is illustrated diagramatically in Fig. 4.

Example A modification of the first-order system given by equations (25) and (26) is used to compare the performances of the estimators described in the last section and the estimator described in this section. Here we assume G, J to be uniformly distributed white noise sequences. Figure 5 shows a typical set of simulation results which demonstrate that the estimator described in this

section gives superior performance when the unknown parameters are varying rapidly with time.

211

Adaptive estimation using parallel processing techniques reset ~(~f ~_,,Oi)l.~( ~1 ~_,.8) before updotin O

~ ~ j

I ;('''-''e''

"("

I r.,.* ; , , , ,

0_,-;(,'," - ' e,

_ _ . d ' . ~ ; . . . . . . I "(; ; - , 4 )

I

1 I

I

-J~ *

I

_ / ,o, o,,/

reset ; , , ,

' %- 2.select ,

k - , Iv

,8,-;,,

*

,_,

, e,I

before updet ng

~

I~:~;'~T-

I ;(tlt o)

; ( t ~_ fie.)

denotes

;((tit

,

I

,-,

,8)

~[elfk ~-I' O(7~k--I)lO'8(ti/r-Z)]

Fig. 4. Discrete-time adaptive estimator for time-varying 0 [assuming scalar measurements and H[t~lO(tk)] = H(t~)l.

0-2 /~Proposed

(Non

resetting)

4-

o.I

Proposed resetting

o

estimator

L

g

2 I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

J

2 3 4 5 6 7 8 9 IO ii 12 13 14 15 16 17 18 192o Time

Fig. 5. Mean square error of state estimates for time-varying noise covariances [G, J random variables uniformly distributed in (0, 1)].

It might be thought that the scheme proposed in this section could be better than that of the previous section for the example of the previous section since the unknown parameter is really time-varying although slowly varying in this example. Simulations showed that when the parameter variations are slowly varying the non-resetting scheme of the previous section works best. IV. C O N T I N O U S - T I M E R E S U L T S

The 'maximum likelihood' technique described in this paper has the advantage that even for continuous-time measurements, an estimator can be obtained which does not require the exact knowledge of the measurement noise covariance.

212

P. K, S. TAM and J. B. MOORE

T i m e - i n v a r i a n t unknown p a r a m e t e r s

Consider the system

i ( t ) = F ( t l t o ) x (t) +

y(t)

G(tlto)u(t)

t_>O

= H ( t l tO) x (t)

(35) (36)

z(t)=y(t)+v(t)

(37)

where {u(t)}, {v (t)} are independent zero mean Gaussian white noise processes with covariance matrices /, R(t) respectively. The initial state x(to) is a Gaussian random vector with mean ~ (tolto, to) and variance P(to]o, to) and is independent of {u(t)} and {v(t)}. The measurement noise covariance R(t) is assumed to be of the form R(t) = aRo(t) where Ro(t) is known but a is an unknown scalar. The entities F(tlto), G(tlto), H(t[to), £(tolto, to), P(totto, tO) are known if the parameter vector to is known, possibly unknown if to is not knowr,. If to and a are specified, the conditional mean state estimate £,(tlt, tO, a) is given from the

:

SO tho~

F

L

i ~

J

I>_,{e,I,} f o r 011

Fig. 6. Continuous-tuneadaptiveestimatorfor time-mvariantO.

I o* t-t,

t !

J~(tl t,ON)

--

Fig. 7. Continuous-timeadaptiveestimator for time-varying0,

213

Adaptive estimation using parallel processing techniques

Kalman-Bucy filter equations:

Yc(tlt, t#, a) = F(tl~):~(tlt, ~, a) + K(tlt~p, a)[z(t) - p(tlt, ~, a)]

(38)

P(tlt, ~, a) = H(tl~)Y,(tlt, ~, a)

(39)

K (tl~,, c,) = P(tlt, ~, cOH'(tl~)[aRo(t)] -I

(40)

P(tlt, ~, a) = F(tl~)P(tlt, ~,, a) + P(tlt, O, a)F'(tl~ ) - aK(tlO, a)Ro(t)K'(tl~, a) + tT(tl~0)G'(tl~).

(41)

If a is specified but ~ is unknown, then assuming ~, the space of admissible values of ~, is discrete or suitably quantized to a finite number of grid points {qJl..... ~p} with known or assumed a priori probability p (~,lt, a) for each ~, the conditional mean state estimate ~ (tit, a) is given by [5]: p

~(tlt, a) = ~, ~(tlt, ¢,,, a)p(ff, lt, a)

(42)

i=l

p (~,lt, a) = c exp {l(6,lt, a)}p(~,lto, a)

(43)

where c is a normalizing constant independent of ~b~and l(~b~It, a) is the log likelihood function of ~ given a and the measurements {z(l")l~'E(to, t)}. Again, allowing for more parallel processing effort, it might appear to the reader that a natural extension of the above results for the case when a is unspecified is 'readily achieved' using the following approach. Assume A, the range of admissible values of a, can be approximated by the finite set {a, ..... aM}. Then, with known or assumed a priori probability p(~,, ajlto) for each (0,, aj), the conditional mean state estimate would be given from P

M

~(tlt) = ~ ~ Y,(tlt, qJ,,a,)p(O,, a, lt).

(45)

i = i 1=]

However, further difficulty arises because there is no direct and simple technique to recursively update p (~, aj[t). We now propose a simple estimation scheme to yield a suboptimal minimum mean square error estimate for the case when a is unknown. Following the approach described in section II we shall now consider directly the space 0 of hypotheses of possible Kalman filter configurations. For each 0 in O, the conditional mean state estimate P,(t[t, O) is given by the following equations:

~c(tlt, O) = F(tlo)~(tlt, o)+ K(tlO)[z(t)- P(tlt, 0)]

(46)

P(tlt, O) = n(tlo)Y~(tJt, o).

(47)

Again, we approximate O with a finite set of grid points {0, . . . . , ON}, and we make the simplifying assumption that each 0, allows the same range of admissible values of a. Then, for any admissible value of a, the log likelihood function l(01]t, a) of 0i given the measurements {z(T)lz~(to, t)} can be expressed as l(0,lt, a) = L(O, lt)/a

(48)

L(0, lt) = f ' ~'(I"IT, 0,)Ro-l(~')Z(¢)de-~ f l IlY(¢t~"00H~-~de}"

(49)

where

Obviously, given the measurements {z0")l¢~(to, t)}, the most likely 0~, denoted 0, is the 0~ which maximizes L (O~lt), which is independent of a. This observation suggests that in the event

214

p. K. S. TAU and J. B. MOORE

that a is not k n o w n , a very simple state estimator is one yielding the estimate g(tlt, 0). Such an estimator is sketched diagramatically in Fig. 5. Simulation results have not been obtained for this estimator since it is not anticipated that they can extablish any more than those for the discrete-time results.

Time-varying unknown parameters Here we consider obtaining estimators when the u n k n o w n parameter vector O(t) and the m e a s u r e m e n t noise intensity a are modelled as piecewise constant functions of time with 0 ( t ) = O(tk), a ( t ) = a(tk) for tk-< t < tk,~. We further assume that O(tk) can take one of the values {01 . . . . . ON}. The estimation scheme that we propose is as follows. At t = t~, 0(to) is found such that L[O(to)[tl] = max {L[O,(to)[t,]}. The states of all the N K a l m a n filters are then reset to ~'[tllt,, O(to)], viz, the state of the K a l m a n filter conditioned on 0(to). In general, at t = t~.~, 0(tk) is found such that L[0(tk)] ..... O(tk-j)] = max {L[O,(tk)lt~÷~, (~(tk -0]}, where

(9(tk 1) = {O(to) ..... O(tk ,)} and

L[O,(t~)lt, O(t~-l)] =

P[~'lr, (~(tk-0, O(tk)lRo-l(~')Z('C) d r - ~ k

I[~(r[r, O(tk-,), 0(t,))llL-, dr. k

Simulations have not been carried out for this scheme since as for the previous not too much additional information can be gained from this. However, it could be stressed that the important point here is that the ideas investigated in this paper, in contrast to the standard Bayesian approach, can in f a c t be applied to continuous-time systems in a straight/orward manner. REFERENCES I. R. E. Kalman, A new approach to linear filtering and prediction problems, Trans. ASME. J. Basic Engng 82, 34-45 (March 1%0). 2. R. E. Kalman and R. S. Bucy, New results in linear filteringand prediction theory, Trans. ASMEZ Basic Engng $3(D) 95-107 (December 1%1). 3. D. T. Magill, Optimal adaptive estimation of Sam led stochastic processes IEEE Trans. Aut. Cont. AC-10,4, 215-218 (October I%5). 4. F. L. Sims, D. G. Lainiotis and D. T. Magill, Recursive algorithm for the calculation of the adaptive Kalman filter weighting coefficients, IEEE Trans. Aut. Cont. (Corresp.) AC-14 (April 1%9). 5. D. G. Lainiotis, Optimal adaptive estimation: structure and parameter adaptation, IEEE Trans. Aut. Cont. AC-16, 2, 160-169 (April 1971). 6. D. L. Alspach and A. Abiri, Optimal nonlinear estimation for a linear system with unknown plant and measurement noise covariances, Proc. o/ 3rd Symposium on Nonlinear Estimation Theory.and its Applications, San Diego, 11-13(September 1972). 7. D. L. Alspach, A parallel processing solution to the adaptive Kalman filteringproblem with vector measurements. Int. J. Comput. Elect Engng. 1, 1, 83-94 (June 1973). 8. D. L. Alspach, A Bayesian nonlinear filtering algorithm for linear systems with unknown time varying noise statistics, Proc. IEEE Conference on Decision and Control incl. 12th Symposium on Adaptive Processes, 533-539(December 1973). 9. C. G. Hilborn, Jr. and D. G. Lainiotis, Optimal adaptive filter realizations for sample stochastic processes with an unknown parameter, IEEE Trans. Aut. Contr. AC-14, 767-770 (December 1969). 10. G. A. Ackerson and K. S. Fu, On state estimation in switching environments, IEEE Trans. Aut. Contr. AC-15, 10--16 (February 1970).