- Email: [email protected]
State estimation for nonlinear discrete dynamic systems with missing observations
State Es tima tion f or Nonlinear Discrete Dynamic Systems with Missing Obseruations t hl,KERiMDEMiRBA$ University Science
of Illinois at Chicago, (M/C
1.541, P.O.
Department
Box 4348, Chicago,
of Electrical IL 60680,
Engineering
and Computer
U.S.A.
ABSTRACT : State estimation of nonlinear dynamic systems is performed by first approximating dynamic systems by Jinite state machines having trellis diagram representations and observation models by functions interpolating aoailable obseraations in the vicinity of the set of times for which observations are missing. The states are then estimated by using the Viterbi decoding algorithm of Information Theory.
I. Introduction Great efforts have been made by researchers on the recursive state estimation of dynamic systems (or models) with observations which are available for all times in an observation interval. As a result, many estimation schemes, referred to as the classical estimation schemes, have been proposed for these dynamic systems, called the classical dynamic systems (or models) (l&6). The proposed estimation schemes have also been used for practical systems (7). A major obstacle with these classical estimation schemes, such as the extended Kalman filter, is that their dynamic models must be linear functions of the disturbance noise and (additive) observation noise (14). Recently, recursive state estimation of more general dynamic models with observations which are available for all times in an observation interval has been considered in (8,9). These dynamic models can be nonlinear functions of the states, observation noise and disturbance noise, whereas the classical estimation models must be linear with respect to these noises. Thus, the proposed schemes are superior to the classica estimation schemes, such as the extended Kalman filter. This paper extends the approach presented in (9) to the state estimation for dynamic models with missing observations. II. Problem Statement Consider
the discrete dynamic
models defined by
x(I+ 1) = f (1,x(l), w(l)) the state model, z(l) = g(Z, x(l), u(o) the observation t This work was carried out while the author
(I) model,
was visiting Bilkent University,
(2) Ankara,
Turkey.
J:, The Franklin
Institute
0016-0032!90
$3.00+0
00
49
Kerim Demirhq where 1 is a discrete moment in time ; z(I) is an m x 1 observation vector at time 1; v(L) is an n x 1 observation noise vector at time I with zero mean and known statistics; II-([) is an s x 1 disturbance noise vector at time 1 with zero mean and known statistics ; x(0) is an r x I random initial state vector with known statistics ; and x(l), 1 >O, is an r x 1 state vector at time 1; J‘(1, x(l), w(2)) and g(1, x(l), v(l)) are either linear or nonlinear functions which define the state at time I+ 1 and observation at time 1 in terms of the disturbance noise, observation noise and state at time 1; and the initial state x(0) and all samples of the disturbance noise and observation noise are independent. This paper treats the state estimation of dynamic models of (1) and (2) with missing observations, that is, the state sequence from time zero to time L, XL g {x(O), _x(l), . . . , x(L)}, is estimated by using available observations from time one to time L. The set of times for which observations are available and the set of times for which observations are missing are denoted by STOA and STOM, respectively. Hence, the union of STOA and STOM is the set of all times from time one to time L, that is TL g (1,2, . , L). Without loss of generality, it is assumed that the first and last observations ~(1) and z(L) are available. In other words, STOM is a subset of the time sequence from time 2 to time L- 1.
III. Estimation
Scheme
State estimation is performed by first approximating the state model of (1) by a finite state model (or machine) and the observation model of (2) by an approximate observation model. This finite state model is represented by a trellis diagram. Then the Viterbi decoding algorithm (10, 11) is used to estimate the state sequence. The state model and observation model are, respectively, approximated by the finite state model (or machine) defined by
and the approximate
observation z(1) A
model defined by
g(1, x,(l), u(l))
if 1~ STOA
8(1)
if IESTOM
(4)
where g(l) is a function which interpolates observations which are available in the vicinity of STOM (12) ; x,,(O) is a discrete random initial state vector with r,, possible values which approximates the initial state x(O), and the possible values of x,(O) are denoted by s,,(O), X&O), . . , and X,,.,,(O), which are called the initial quantization levels (or the quantization levels at time zero); and s,(l), 1 > 0, is the quantized state at time 1, the number of possible values of which is (say) r,, and the possible values of x,(l) are denoted by s,,(l), _Y,,~(]),. , and .~~,,,(l). In the finite state model of (3) : IV,,(~)is a discrete disturbance noise vector with s, possible values which approximates the disturbance noise vector at time 1, and the possible values of ,r,,(l) are denoted by or,,,, ~‘,,~(l),. , and N!,,,,(L); and Q( .) is the quantizer defined in (8). This quantizer is a function which first divides the r-dimensional Euclidean space into nonoverlapping generalized rectangles of the same size (called
State Estimation.for Time zero
Time 1
Time i-l
Nonlinear Dynamic Systems
Time
I
Time L
FIG. 1. Trellis diagram of state. the gate size) such that the union of all rectangles is the r-dimensional Euclidean space. It then assigns to each rectangle (called the gate) the center point of the rectangle. The finite state model of (3) is represented by a trellis diagram, Fig. 1, from time zero to time L (called the trellis diagram of the state), and the following metrics are assigned to the nodes, branches and paths of this trellis diagram. The metric of the branch connecting a node (or quantization level) xy,(Z- 1) to another node x,,(l) is denoted by MB(xy,(lI) +x&l)) and defined by MWx,, (l-
1) + &pi (0) A In {Wd-
1) + ~,,(0>~(40lx,,(0~1
where In indicates the natural logarithm; p(z(l)Ix,,,(l)) is the conditional density function of the observation r(l) given that x(l) = x,,,,!(o (in Eq. 2), where if Z(I) is missing, it is first estimated by j(I) and then used in this conditional density probability from the node function ; and T(x,,,(l- 1) + x ‘,,,1(r)) is the transition x,,,,(I- 1) to the node x&l>, which is easily determined by the finite state model of (3). The metric of a node x4,,J1) is denoted by MN(x,,,,(l)) and defined by MN(x,,,,(r))
g
In {Bob o
{x,(r) = x,~(/)>}
if I = 0 otherwise
The metric of a path is defined as the sum of all metrics of the nodes and branches along the path. The state estimation problem is to find a path through the trellis diagram along which the quantization levels will be the estimate of the state sequence. This path is found by using multiple hypothesis testing (2). It can be shown (8) that the optimum testing rule, which minimizes the overall error probability per decision, is to choose the path with the greatest metric (and in the case that there are more
Kerim DemirbaJ
than one path having the same greatest metric, then to choose anyone of these paths having the same greatest metric). The path with the greatest metric is found by the Viterbi decoding algorithm (VDA).
IV. Performance The performance of the proposed scheme, which is based upon the performance of the VDA used to estimate the states, can be quantified by a Gallager-type ensemble upper bound (11, 13), since the evaluation of the exact error probability or error probability bound for choosing the correct path is complex. It can be shown (8) that such an ensemble bound is given by
P, < G(v) fj /= I
c d4PMl)l4’~“”
XEX’
1 i I+
Y
dz(l)
for any VE [0, 11,
(5)
with G(v) A (K-l)@($)‘:‘+“] where K is the number of possible paths through the trellis diagram of the state; 7~;“’ and rcy are the minimum and maximum values of the occurrence probabilities of possible values of the discrete initial state x,(O) ; $‘I” and rc;lax, I > 0, are the minimum and maximum values of the transition probabilities from time I- 1 to time 1 (in the trellis diagram of the state), respectively; F’ is the set of all possible quantization levels of the states from time one to time L; P, is the ensemble averaged overall error probability for state estimation ; q( .) is an arbitrary probability density function on _Y ; and p(z(Z) Ix) is the conditional density function of the observation at time 1 given that x(l) = x. As the performance measure of the proposed scheme, the uniformly weighted ensemble bound with v = 1 in (5) is used since it is the easiest bound to evaluate weighted” means that q(x) = l/M’ for all x, where M’ is the (8) ; “ uniformly number of elements in _Y’.Consider, as an example, the models defined by x(l+
1) = f(l, x(l), w(l)) the state model z(l) = g(l, x(o) +v(l)
the observation
(6) model
(7)
where x(0) and v(l) are assumed to be Gaussian noises with means m,,, 0; and covariances R,, and R,(I), respectively. Substituting p(z(l)]x(l) = x), v = 1, and q( .) = (l/M’) into the bound in (5), the following bound is obtained (8) :
where
W,xl, 52
x2) A exp I- MLxd -dLx2>lTK ‘(I)MLxd -_g(Axdl/~~ Journalof
the Franklin Pergamon
lnstltute Press plc
State Estimation for Nonlinear
Dynamic
Systems
where the superscript T indicates the transpose. The bound of (8) is the one used as the performance measure of the proposed scheme for the models of (6) and (7).
V. Simulations State estimation of many dynamic systems with white Gaussian noise and missing observations were simulated on the IBM 3081K mainframe computer. The disturbance noise and initial state were approximated by the discrete random variables given in (8). The approximating discrete random variables were also assumed to be time-invariant. Simulation results of three examples are presented in Figs 2(a)4(c). In these examples. it was assumed that the set of available observations (SAO) is {z(l), z(2), z(5), z(6)} and the set of missing observations (SMO) was {z(3), z(4)}. The state sequence {x(O), X(I), x(2), x(3), x(3), x(4), x(5), ~(6)) were estimated by using SAO. This estimation was performed by both the proposed scheme and (extended) Kalman filter. The interpolating function g(r) in (4) was taken to be a polynomial which interpolates available observations in the vicinity of STOM (12). In Figs 2(a)-4(c) the simulated state and observation models are given at the first and second rows at the top left-hand corner; ACTUAL, KALMAN and OPDS represent the actual and estimated values of the states by the (extended) Kalman filter and the proposed scheme ; AAEK and AAEOP denote the time XIK+l1=1.SX[KI+WIKI ZIKl:ZXIK)+VIKI NUil. OF DISC. vRRlXI011:1.000
0
zz &B t--Lo
FOR
X[O1=3
NUN. OF DISC. FOR VRRlW[.11=2.700 VRRlVl.lI:O.SOO GATE SlZE:0.25C
W[.l=3
LEGEND 0:
RCTURL
A:
KALtlRN
+:
OPOS
FIG. 2. (a) Actual and estimated values of states. Vol.327,No. ,.,q.49-59.1990 Prmted in Great Britain
53
X,K+,I=l.5X[KI+WIK~ ZIKl=ZX[Kl+VIKl NUfl. OF DISCVRRLXIOI~~l.ODO
FOR
ElXlDll-1.300 NUtI. OF DISC. FOR VRR[WI.11=2.700 VRRlVl.1lz0.600 GATE 5IZE=0.250
X101=3
LEGEND
Wl.1=3
0:
ER.COV.
EDUNO=0.5144ZE-3
00
1.20
2.40
FIG. 2. (b) Error variances
TIME
3.60
/
I
4.80
6.00
and bound for estimates
of states.
XIK+11=1.5XlKltWIK1 2IKI~PXlKI+VIK,
z i, L
I
NUM. OF OISC- FOR VRRLXIOI~=I.OOO EIXlOIl=1.300 NUM. OF DISC. FOR VARIWI .11=2.700 VRR~VI.l~=D.BOO GRTE
X101=3
S12E=0.250
=:
54
A:
KRLnAN
t:
OPOS
RREK=D.B6297ZEO RREOP=O.IOB785EI
cc
FIG. 2. (c) Absolute
LEGEND
wL.1~3
and time averaged
absolute
errors for estimates
of states.
Journal of the Frankhn Institute Pergamon Press plc
XIK+lI:X~KlEXPIO.OZXlKII+W~Kl ZlKI:4X~KI+VIKI NUH. DF DISC. FOR VRR[X~01)=1.000 ElX~011=2.000 NUFI. OF DISC. FOR VRR[WI. ll=Z.OOO VRR~VI.Il:l.OOO GATE
XlO1=3
LEGEND
W[.lz3
SIZE=O.ZOO
1.20
00
3.60
2.40
0:
RCTURL
A:
KALIIRN
+:
OPOS
6.00
4.80
TIME
FIG. 3. (a) Actual
and estimated
NUM. OF DISCVRRlX1011=1.000 EIX~OIl=Z.OOO NUN. OF DISC. VRR~WI.11=2.0DO VRRlV1.11:1.000 GRTE
FOR
X101=3
FOR
WC.123
values of states.
LEGEND c]:
ER.COV.
SlZE=O.ZOO BDUNO=O.l5996E-1
t$
1.20
FIG. 3. (b) Error variances Vol 327,No. I,pp.49 59, 1990 Printed in Great Braan
2.40
TIME
3.60
4.80
and bound for estimates
6.00
of states.
55
XIK+lI=XIKlEXPIO.DPXIKIl+WIK1 ZIK1~4X~KlrVIKl NUil. OF DISC. FDR VRR~XIOl)=l.OOO EIX1011=2.000 NUM. OF DISC. FOR VRRINI.11=2.000 VRR~VI.11=1.000 GRTE 51ZE=0.200
X101=3
Wc.1~3 A:
KRLnRN
+:
OPOS
RREK=O.S39526EO RREOP=O.513728EO
so
FIG. 3. (c) Absolute
1
.zo
2.40
TIME
and time averaged
3.60
absolute
4.80
6.00
errors for estimates
of states
X~K+LI:I.ZXIKI+WIKI Z[Kl=O.EXfKI+VlKl NUtI. OF CIISC. FOR VAR1XI011~I.000 EIXl011:4.000 NUM. OF DISC. FOR VRRIWI.ll=3.000 VRA[YI.1)~0.500 GRTE SIZE=O.ZOO
I
XlO1=3
LEGEND
Wt.1:3
0:
RCTURL
A:
KRLflRN
+:
OPOS
: ; :
$.00
1.20
2.40
TIME
3.60
FIG. 4. (a) Actual and estimated
4.80
1
6.00
values of states. Journal of the Franklm lnst~tute Pergamon Press plc
XlK+il=l.ZXlKI+WlKI ZIKI:O.BX~Kl+VIKI NUM. DF OISC- FOR VRR~XIOII=l.OOO EIX~OlI=4.ODO NUN. OF DISC. FOR VAR[CIf .11=3.000 VRRIVl.ll~0.500 GATE
X101=3
LEGEND
W[.l=3
q:
ER-COV
SIZE=O.ZOO BOUNO=0.63566E-4
TINE
FIG. 4. (b) Error variances
XlK+l1:1.EX~KI+WIKI ZIKI=O.BXfKl+VIKl NUFI. OF DISC. FOR !!RRlX~01):1.000 ElXIOll=4.000 NUM. OF OISC. FOR VARlUI.11=3.000 VRRlVL.1~=0.500 GATE
and bound for estimates
of states.
XcO1.3
LEGEND
Wc.1~3
SlZE=O.ZOO
A:
KRL,'lRN
+:
OPOS
RAEK:O.153433E1 RREOP=O.S231lOEO
FIG. 4. (c) Absolute Vol.327.No. l.pp.49%59,1990 Printedin Great Elnta~n
and time averaged
absolute
errors for estimates
of states.
57
Kerim
Demirha?
averaged absolute errors for the estimates obtained by the (extended) Kalman filter and proposed scheme; E(A) and VAR(A) indicate the mean value and variance of the random variable A, respectively ; GATE SIZE shows the gate size used for the quantizer in (3) ; BOUND represents the bound of (8) which was used as the performance measure of the proposed scheme; ER. COV. stands for the error variances for the estimates obtained by the (extended) Kalman filter; and NUM. OF DISC. FOR A shows the number of possible values of the discrete random variable which was used to approximate the random variable A. For Figs 2(a)-3(c) : the interpolating function @(r>in (4) was taken to be a firstorder polynomial which interpolates the observations at times 2 and 5, whereas J(1) was taken to be a second-order polynomial which interpolates the observations at times 2, 5 and 6 for Fig. 4(a)-(c). The estimates obtained by the proposed scheme are better than the Kalman estimates for the nonlinear dynamic models given in Figs 3(a)4(c). On the other hand, the Kalman estimates are better for the linear dynamic mode1 in Fig. 2(a))(c), since the Kalman filter is optimum for linear dynamic models with white Gaussian noise. The bound of (8) may sometimes become a number greater than one (i.e. unless) for some dynamic models since it is derived by using several inequalities. One also must keep in mind that this is an ensemble bound for the performance of the proposed scheme, which does not exactly yield the performance of the proposed scheme for a specific dynamic system (8).
VI. Conclusions An estimation scheme is presented for nonlinear dynamic systems with missing observations. This estimation scheme uses interpolating functions. The performance of this scheme is based upon the approximation accuracy of the interpolating functions used for missing observations. The proposed estimation scheme is superior to the classical estimation schemes, such as the extended Kalman filter, since the models of the proposed scheme (which can be nonlinear functions of the states, disturbance noise, and observation noise) are more genera1 than the models of the classical estimation schemes (which must be linear functions of the disturbance noise and observation noise). References (1) R. E. Kalman, “A new approach to linear filtering and prediction problems”, Trans. ASME J. Basic Engng, ser. D., Vol. 82, pp. 35545, 1960. (2) A. P. Sage and J. L. Melsa, “Estimation Theory with Applications to Communications and Control”, McGraw-Hill, New York, 1971. (3) T. Kailath, “An innovation approach to least square estimation, part 1 : linear filtering in adaptive white noise”, IEEE Trans. Aut. Control, Vol. AC-13, No. 6, pp. 646 655, 1968. (4) T. Kailath, “A view of three decades in linear filtering theory”, IEEE Trans. Inj: Theory, Vol. IT-20, No. 2, pp. 146-181, 1974. (5) J. Makhoul, “Linear prediction : a tutorial review”, Proc. IEEE, Vol. 63, pp. 561-580, 1975.
58
Journal
of the Franklm Pergamon
Institute Press plc
State
Estimation,for
Nonlinear
Dynamic
Systems
(6) J. S. Medich, “A survey of data smoothing for linear and nonlinear dynamic systems”, Automatica, Vol. 9, pp. 151-162, 1973. (7) C. E. Hutchinson, “The Kalman Filter applied to aerospace and electronic systems”, IEEE Trans. Aero. Electron. Systems, Vol. AES-20, No. 4, pp. 500-504, 1984. (8) K. Demirbas, “Information theoretic smoothing algorithms for dynamic systems with or without interference”, in Advances in Control and Dynamic Systems, Volume XXI, Academic Press, New York, pp. 175-295, 1984. (9) K. Demirbas and C. T. Leondes, “Optimum decoding based smoothing algorithm for dynamic systems”, ht. J. Systems Sci., Vol. 16, No. 8, pp. 951-966, 1985. (10) G. D. Forney, Jr, “Convolutional Codes II. Maximum Likelihood Decoding”. It$ Control, Vol. 25, pp. 222-266, 1974. (11) A. J. Viterbi and J. K. Omura, “Principles of Digital Communication and Coding”, McGraw-Hill, New York, 1979. (12) S. D. Conte and Carl de Boor, “Elementary Numerical Analysis”, McGraw-Hill, New York, 1972. (13) R. G. Gallager, “A simple derivation of the coding theorem and some applications”, IEEE Trans. If. Theory, Vol. IT-1 1, pp. 3-18, 1965.
Vol. 327, No. I, pp. 49-59, 1990 Prmtcd m Great Britain
59
of Illinois at Chicago, (M/C
1.541, P.O.
Department
Box 4348, Chicago,
of Electrical IL 60680,
Engineering
and Computer
U.S.A.
ABSTRACT : State estimation of nonlinear dynamic systems is performed by first approximating dynamic systems by Jinite state machines having trellis diagram representations and observation models by functions interpolating aoailable obseraations in the vicinity of the set of times for which observations are missing. The states are then estimated by using the Viterbi decoding algorithm of Information Theory.
I. Introduction Great efforts have been made by researchers on the recursive state estimation of dynamic systems (or models) with observations which are available for all times in an observation interval. As a result, many estimation schemes, referred to as the classical estimation schemes, have been proposed for these dynamic systems, called the classical dynamic systems (or models) (l&6). The proposed estimation schemes have also been used for practical systems (7). A major obstacle with these classical estimation schemes, such as the extended Kalman filter, is that their dynamic models must be linear functions of the disturbance noise and (additive) observation noise (14). Recently, recursive state estimation of more general dynamic models with observations which are available for all times in an observation interval has been considered in (8,9). These dynamic models can be nonlinear functions of the states, observation noise and disturbance noise, whereas the classical estimation models must be linear with respect to these noises. Thus, the proposed schemes are superior to the classica estimation schemes, such as the extended Kalman filter. This paper extends the approach presented in (9) to the state estimation for dynamic models with missing observations. II. Problem Statement Consider
the discrete dynamic
models defined by
x(I+ 1) = f (1,x(l), w(l)) the state model, z(l) = g(Z, x(l), u(o) the observation t This work was carried out while the author
(I) model,
was visiting Bilkent University,
(2) Ankara,
Turkey.
J:, The Franklin
Institute
0016-0032!90
$3.00+0
00
49
Kerim Demirhq where 1 is a discrete moment in time ; z(I) is an m x 1 observation vector at time 1; v(L) is an n x 1 observation noise vector at time I with zero mean and known statistics; II-([) is an s x 1 disturbance noise vector at time 1 with zero mean and known statistics ; x(0) is an r x I random initial state vector with known statistics ; and x(l), 1 >O, is an r x 1 state vector at time 1; J‘(1, x(l), w(2)) and g(1, x(l), v(l)) are either linear or nonlinear functions which define the state at time I+ 1 and observation at time 1 in terms of the disturbance noise, observation noise and state at time 1; and the initial state x(0) and all samples of the disturbance noise and observation noise are independent. This paper treats the state estimation of dynamic models of (1) and (2) with missing observations, that is, the state sequence from time zero to time L, XL g {x(O), _x(l), . . . , x(L)}, is estimated by using available observations from time one to time L. The set of times for which observations are available and the set of times for which observations are missing are denoted by STOA and STOM, respectively. Hence, the union of STOA and STOM is the set of all times from time one to time L, that is TL g (1,2, . , L). Without loss of generality, it is assumed that the first and last observations ~(1) and z(L) are available. In other words, STOM is a subset of the time sequence from time 2 to time L- 1.
III. Estimation
Scheme
State estimation is performed by first approximating the state model of (1) by a finite state model (or machine) and the observation model of (2) by an approximate observation model. This finite state model is represented by a trellis diagram. Then the Viterbi decoding algorithm (10, 11) is used to estimate the state sequence. The state model and observation model are, respectively, approximated by the finite state model (or machine) defined by
and the approximate
observation z(1) A
model defined by
g(1, x,(l), u(l))
if 1~ STOA
8(1)
if IESTOM
(4)
where g(l) is a function which interpolates observations which are available in the vicinity of STOM (12) ; x,,(O) is a discrete random initial state vector with r,, possible values which approximates the initial state x(O), and the possible values of x,(O) are denoted by s,,(O), X&O), . . , and X,,.,,(O), which are called the initial quantization levels (or the quantization levels at time zero); and s,(l), 1 > 0, is the quantized state at time 1, the number of possible values of which is (say) r,, and the possible values of x,(l) are denoted by s,,(l), _Y,,~(]),. , and .~~,,,(l). In the finite state model of (3) : IV,,(~)is a discrete disturbance noise vector with s, possible values which approximates the disturbance noise vector at time 1, and the possible values of ,r,,(l) are denoted by or,,,, ~‘,,~(l),. , and N!,,,,(L); and Q( .) is the quantizer defined in (8). This quantizer is a function which first divides the r-dimensional Euclidean space into nonoverlapping generalized rectangles of the same size (called
State Estimation.for Time zero
Time 1
Time i-l
Nonlinear Dynamic Systems
Time
I
Time L
FIG. 1. Trellis diagram of state. the gate size) such that the union of all rectangles is the r-dimensional Euclidean space. It then assigns to each rectangle (called the gate) the center point of the rectangle. The finite state model of (3) is represented by a trellis diagram, Fig. 1, from time zero to time L (called the trellis diagram of the state), and the following metrics are assigned to the nodes, branches and paths of this trellis diagram. The metric of the branch connecting a node (or quantization level) xy,(Z- 1) to another node x,,(l) is denoted by MB(xy,(lI) +x&l)) and defined by MWx,, (l-
1) + &pi (0) A In {Wd-
1) + ~,,(0>~(40lx,,(0~1
where In indicates the natural logarithm; p(z(l)Ix,,,(l)) is the conditional density function of the observation r(l) given that x(l) = x,,,,!(o (in Eq. 2), where if Z(I) is missing, it is first estimated by j(I) and then used in this conditional density probability from the node function ; and T(x,,,(l- 1) + x ‘,,,1(r)) is the transition x,,,,(I- 1) to the node x&l>, which is easily determined by the finite state model of (3). The metric of a node x4,,J1) is denoted by MN(x,,,,(l)) and defined by MN(x,,,,(r))
g
In {Bob o
{x,(r) = x,~(/)>}
if I = 0 otherwise
The metric of a path is defined as the sum of all metrics of the nodes and branches along the path. The state estimation problem is to find a path through the trellis diagram along which the quantization levels will be the estimate of the state sequence. This path is found by using multiple hypothesis testing (2). It can be shown (8) that the optimum testing rule, which minimizes the overall error probability per decision, is to choose the path with the greatest metric (and in the case that there are more
Kerim DemirbaJ
than one path having the same greatest metric, then to choose anyone of these paths having the same greatest metric). The path with the greatest metric is found by the Viterbi decoding algorithm (VDA).
IV. Performance The performance of the proposed scheme, which is based upon the performance of the VDA used to estimate the states, can be quantified by a Gallager-type ensemble upper bound (11, 13), since the evaluation of the exact error probability or error probability bound for choosing the correct path is complex. It can be shown (8) that such an ensemble bound is given by
P, < G(v) fj /= I
c d4PMl)l4’~“”
XEX’
1 i I+
Y
dz(l)
for any VE [0, 11,
(5)
with G(v) A (K-l)@($)‘:‘+“] where K is the number of possible paths through the trellis diagram of the state; 7~;“’ and rcy are the minimum and maximum values of the occurrence probabilities of possible values of the discrete initial state x,(O) ; $‘I” and rc;lax, I > 0, are the minimum and maximum values of the transition probabilities from time I- 1 to time 1 (in the trellis diagram of the state), respectively; F’ is the set of all possible quantization levels of the states from time one to time L; P, is the ensemble averaged overall error probability for state estimation ; q( .) is an arbitrary probability density function on _Y ; and p(z(Z) Ix) is the conditional density function of the observation at time 1 given that x(l) = x. As the performance measure of the proposed scheme, the uniformly weighted ensemble bound with v = 1 in (5) is used since it is the easiest bound to evaluate weighted” means that q(x) = l/M’ for all x, where M’ is the (8) ; “ uniformly number of elements in _Y’.Consider, as an example, the models defined by x(l+
1) = f(l, x(l), w(l)) the state model z(l) = g(l, x(o) +v(l)
the observation
(6) model
(7)
where x(0) and v(l) are assumed to be Gaussian noises with means m,,, 0; and covariances R,, and R,(I), respectively. Substituting p(z(l)]x(l) = x), v = 1, and q( .) = (l/M’) into the bound in (5), the following bound is obtained (8) :
where
W,xl, 52
x2) A exp I- MLxd -dLx2>lTK ‘(I)MLxd -_g(Axdl/~~ Journalof
the Franklin Pergamon
lnstltute Press plc
State Estimation for Nonlinear
Dynamic
Systems
where the superscript T indicates the transpose. The bound of (8) is the one used as the performance measure of the proposed scheme for the models of (6) and (7).
V. Simulations State estimation of many dynamic systems with white Gaussian noise and missing observations were simulated on the IBM 3081K mainframe computer. The disturbance noise and initial state were approximated by the discrete random variables given in (8). The approximating discrete random variables were also assumed to be time-invariant. Simulation results of three examples are presented in Figs 2(a)4(c). In these examples. it was assumed that the set of available observations (SAO) is {z(l), z(2), z(5), z(6)} and the set of missing observations (SMO) was {z(3), z(4)}. The state sequence {x(O), X(I), x(2), x(3), x(3), x(4), x(5), ~(6)) were estimated by using SAO. This estimation was performed by both the proposed scheme and (extended) Kalman filter. The interpolating function g(r) in (4) was taken to be a polynomial which interpolates available observations in the vicinity of STOM (12). In Figs 2(a)-4(c) the simulated state and observation models are given at the first and second rows at the top left-hand corner; ACTUAL, KALMAN and OPDS represent the actual and estimated values of the states by the (extended) Kalman filter and the proposed scheme ; AAEK and AAEOP denote the time XIK+l1=1.SX[KI+WIKI ZIKl:ZXIK)+VIKI NUil. OF DISC. vRRlXI011:1.000
0
zz &B t--Lo
FOR
X[O1=3
NUN. OF DISC. FOR VRRlW[.11=2.700 VRRlVl.lI:O.SOO GATE SlZE:0.25C
W[.l=3
LEGEND 0:
RCTURL
A:
KALtlRN
+:
OPOS
FIG. 2. (a) Actual and estimated values of states. Vol.327,No. ,.,q.49-59.1990 Prmted in Great Britain
53
X,K+,I=l.5X[KI+WIK~ ZIKl=ZX[Kl+VIKl NUfl. OF DISCVRRLXIOI~~l.ODO
FOR
ElXlDll-1.300 NUtI. OF DISC. FOR VRR[WI.11=2.700 VRRlVl.1lz0.600 GATE 5IZE=0.250
X101=3
LEGEND
Wl.1=3
0:
ER.COV.
EDUNO=0.5144ZE-3
00
1.20
2.40
FIG. 2. (b) Error variances
TIME
3.60
/
I
4.80
6.00
and bound for estimates
of states.
XIK+11=1.5XlKltWIK1 2IKI~PXlKI+VIK,
z i, L
I
NUM. OF OISC- FOR VRRLXIOI~=I.OOO EIXlOIl=1.300 NUM. OF DISC. FOR VARIWI .11=2.700 VRR~VI.l~=D.BOO GRTE
X101=3
S12E=0.250
=:
54
A:
KRLnAN
t:
OPOS
RREK=D.B6297ZEO RREOP=O.IOB785EI
cc
FIG. 2. (c) Absolute
LEGEND
wL.1~3
and time averaged
absolute
errors for estimates
of states.
Journal of the Frankhn Institute Pergamon Press plc
XIK+lI:X~KlEXPIO.OZXlKII+W~Kl ZlKI:4X~KI+VIKI NUH. DF DISC. FOR VRR[X~01)=1.000 ElX~011=2.000 NUFI. OF DISC. FOR VRR[WI. ll=Z.OOO VRR~VI.Il:l.OOO GATE
XlO1=3
LEGEND
W[.lz3
SIZE=O.ZOO
1.20
00
3.60
2.40
0:
RCTURL
A:
KALIIRN
+:
OPOS
6.00
4.80
TIME
FIG. 3. (a) Actual
and estimated
NUM. OF DISCVRRlX1011=1.000 EIX~OIl=Z.OOO NUN. OF DISC. VRR~WI.11=2.0DO VRRlV1.11:1.000 GRTE
FOR
X101=3
FOR
WC.123
values of states.
LEGEND c]:
ER.COV.
SlZE=O.ZOO BDUNO=O.l5996E-1
t$
1.20
FIG. 3. (b) Error variances Vol 327,No. I,pp.49 59, 1990 Printed in Great Braan
2.40
TIME
3.60
4.80
and bound for estimates
6.00
of states.
55
XIK+lI=XIKlEXPIO.DPXIKIl+WIK1 ZIK1~4X~KlrVIKl NUil. OF DISC. FDR VRR~XIOl)=l.OOO EIX1011=2.000 NUM. OF DISC. FOR VRRINI.11=2.000 VRR~VI.11=1.000 GRTE 51ZE=0.200
X101=3
Wc.1~3 A:
KRLnRN
+:
OPOS
RREK=O.S39526EO RREOP=O.513728EO
so
FIG. 3. (c) Absolute
1
.zo
2.40
TIME
and time averaged
3.60
absolute
4.80
6.00
errors for estimates
of states
X~K+LI:I.ZXIKI+WIKI Z[Kl=O.EXfKI+VlKl NUtI. OF CIISC. FOR VAR1XI011~I.000 EIXl011:4.000 NUM. OF DISC. FOR VRRIWI.ll=3.000 VRA[YI.1)~0.500 GRTE SIZE=O.ZOO
I
XlO1=3
LEGEND
Wt.1:3
0:
RCTURL
A:
KRLflRN
+:
OPOS
: ; :
$.00
1.20
2.40
TIME
3.60
FIG. 4. (a) Actual and estimated
4.80
1
6.00
values of states. Journal of the Franklm lnst~tute Pergamon Press plc
XlK+il=l.ZXlKI+WlKI ZIKI:O.BX~Kl+VIKI NUM. DF OISC- FOR VRR~XIOII=l.OOO EIX~OlI=4.ODO NUN. OF DISC. FOR VAR[CIf .11=3.000 VRRIVl.ll~0.500 GATE
X101=3
LEGEND
W[.l=3
q:
ER-COV
SIZE=O.ZOO BOUNO=0.63566E-4
TINE
FIG. 4. (b) Error variances
XlK+l1:1.EX~KI+WIKI ZIKI=O.BXfKl+VIKl NUFI. OF DISC. FOR !!RRlX~01):1.000 ElXIOll=4.000 NUM. OF OISC. FOR VARlUI.11=3.000 VRRlVL.1~=0.500 GATE
and bound for estimates
of states.
XcO1.3
LEGEND
Wc.1~3
SlZE=O.ZOO
A:
KRL,'lRN
+:
OPOS
RAEK:O.153433E1 RREOP=O.S231lOEO
FIG. 4. (c) Absolute Vol.327.No. l.pp.49%59,1990 Printedin Great Elnta~n
and time averaged
absolute
errors for estimates
of states.
57
Kerim
Demirha?
averaged absolute errors for the estimates obtained by the (extended) Kalman filter and proposed scheme; E(A) and VAR(A) indicate the mean value and variance of the random variable A, respectively ; GATE SIZE shows the gate size used for the quantizer in (3) ; BOUND represents the bound of (8) which was used as the performance measure of the proposed scheme; ER. COV. stands for the error variances for the estimates obtained by the (extended) Kalman filter; and NUM. OF DISC. FOR A shows the number of possible values of the discrete random variable which was used to approximate the random variable A. For Figs 2(a)-3(c) : the interpolating function @(r>in (4) was taken to be a firstorder polynomial which interpolates the observations at times 2 and 5, whereas J(1) was taken to be a second-order polynomial which interpolates the observations at times 2, 5 and 6 for Fig. 4(a)-(c). The estimates obtained by the proposed scheme are better than the Kalman estimates for the nonlinear dynamic models given in Figs 3(a)4(c). On the other hand, the Kalman estimates are better for the linear dynamic mode1 in Fig. 2(a))(c), since the Kalman filter is optimum for linear dynamic models with white Gaussian noise. The bound of (8) may sometimes become a number greater than one (i.e. unless) for some dynamic models since it is derived by using several inequalities. One also must keep in mind that this is an ensemble bound for the performance of the proposed scheme, which does not exactly yield the performance of the proposed scheme for a specific dynamic system (8).
VI. Conclusions An estimation scheme is presented for nonlinear dynamic systems with missing observations. This estimation scheme uses interpolating functions. The performance of this scheme is based upon the approximation accuracy of the interpolating functions used for missing observations. The proposed estimation scheme is superior to the classical estimation schemes, such as the extended Kalman filter, since the models of the proposed scheme (which can be nonlinear functions of the states, disturbance noise, and observation noise) are more genera1 than the models of the classical estimation schemes (which must be linear functions of the disturbance noise and observation noise). References (1) R. E. Kalman, “A new approach to linear filtering and prediction problems”, Trans. ASME J. Basic Engng, ser. D., Vol. 82, pp. 35545, 1960. (2) A. P. Sage and J. L. Melsa, “Estimation Theory with Applications to Communications and Control”, McGraw-Hill, New York, 1971. (3) T. Kailath, “An innovation approach to least square estimation, part 1 : linear filtering in adaptive white noise”, IEEE Trans. Aut. Control, Vol. AC-13, No. 6, pp. 646 655, 1968. (4) T. Kailath, “A view of three decades in linear filtering theory”, IEEE Trans. Inj: Theory, Vol. IT-20, No. 2, pp. 146-181, 1974. (5) J. Makhoul, “Linear prediction : a tutorial review”, Proc. IEEE, Vol. 63, pp. 561-580, 1975.
58
Journal
of the Franklm Pergamon
Institute Press plc
State
Estimation,for
Nonlinear
Dynamic
Systems
(6) J. S. Medich, “A survey of data smoothing for linear and nonlinear dynamic systems”, Automatica, Vol. 9, pp. 151-162, 1973. (7) C. E. Hutchinson, “The Kalman Filter applied to aerospace and electronic systems”, IEEE Trans. Aero. Electron. Systems, Vol. AES-20, No. 4, pp. 500-504, 1984. (8) K. Demirbas, “Information theoretic smoothing algorithms for dynamic systems with or without interference”, in Advances in Control and Dynamic Systems, Volume XXI, Academic Press, New York, pp. 175-295, 1984. (9) K. Demirbas and C. T. Leondes, “Optimum decoding based smoothing algorithm for dynamic systems”, ht. J. Systems Sci., Vol. 16, No. 8, pp. 951-966, 1985. (10) G. D. Forney, Jr, “Convolutional Codes II. Maximum Likelihood Decoding”. It$ Control, Vol. 25, pp. 222-266, 1974. (11) A. J. Viterbi and J. K. Omura, “Principles of Digital Communication and Coding”, McGraw-Hill, New York, 1979. (12) S. D. Conte and Carl de Boor, “Elementary Numerical Analysis”, McGraw-Hill, New York, 1972. (13) R. G. Gallager, “A simple derivation of the coding theorem and some applications”, IEEE Trans. If. Theory, Vol. IT-1 1, pp. 3-18, 1965.
Vol. 327, No. I, pp. 49-59, 1990 Prmtcd m Great Britain
59