LlNEAR AND NONLINEAR DESCRETE FiLTERiNG
FOR CONTiNUOUS SYSTEMS

Thesis for the Degree of Ph, D.
MICH§GAN STATE UNNERSWY
MARVIN A. NEEDLER
:97;

 

fHESl9 .. LIBRAP .

, -'Michigan State
University

 

This is to certify that the

thesis entitled

LINEAR AND NONLINEAR DISCRETE FILTERING
FOR CONTINUOUS SYSTEMS

presented by

Marvin A . Needler

has been accepted towards fulfillment
of the requirements for

X4244? f/gf/ I

Major professor

Date February 12, 1971

 

O~169

ABSTRACT

LINEAR AND NONLINEAR DISCRETE FILTERING FOR CONTINUOUS SYSTEMS

by
Marvin A. Needler

The subject of this thesis is the synthesis of discrete-time filters
for real-time estimation of the state of continuous-time dynamic systems.

The dynamic system is modeled by a set of stochastic differential
equations and the stochastic disturbance is modeled as a Brownian motion
process. The filtering process is assumed to utilize en-line digital
computation and thus a discrete-time recursiVe filter is developed to
estimate the state of the system based on sampled-data measurements.

Starting with the equations for propagation of the probability den-
sity function of the state conditioned on previous measurements, the ’
filtering equations are developed both for linear and nonlinear filtering.

In the linear filtering case, two new methods are suggested for
reducing computation and increasing accuracy. These methods are illus—
trated in the digital simulation of state estimation of a steam turbo-
generator system. In addition, the methods of divergence compensation
are reviewed.

In the nonlinear filtering case, a second—order filter is suggested,

and it is used to estimate the state of a reentry body. An original

Marvin A. Needler

adaptive filter is mechanized to reduce filter error and is shown to
improve the estimation accuracy of the linear and nonlinear filter in
reentry estimation. The adaptive filter adapts the filter gain to the
statistics of the measurement estimation errors without requiring
significant computation and.without altering the state covariance com-
putation. Based on the reentry estimation problem, the adaptive filter
produces an improved state estimate that approaches the optimum estimate

for the given filtering problem.

LINEAR.AND NONLINEAR DISCRETE FILTERING FOR CONTINUOUS SYSTEMS

By

Era:

Marvin A: Needler

A THESIS

Submitted to
‘Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering and Systems Science

1971

ACKNOWLEDGEMENTS

I would like to express my gratitude to my major professor, Dr.
Gerald L. Park, for the guidance and inspiration that he provided
throughout my doctoral program. In addition, I am deeply grateful to
the other guidance committee members, Dr. Richard C. Dubes, Dr..J.
Sutherland Frame, Dr. William L. Kilmer, and Dr. John B. Kreer, for
their invaluable contributions to my doctoral program. .A special note
should be made of the joint effort provided by Dr. Frame in developing
Section 3-2 of this thesis.

Support was granted by Consumers' Power Company Cooperative Research
Project athichigan State University for one summer and by the National
Science Foundation for one year. This support is gratefully acknowledged.

In addition, I would like to thank Mr. Howard Bell, who introduced
me to filtering theory. Also, I am indebted to Mrs. Ann Darrenkamp for
typing the final manuscript. And finally, the original manuscript was

typed by my wife, Priscilla, whose support made this thesis possible.

ii

TABLE OF CONTENTS

Chapter I Introduction

Chapter II Analysis and Development of Filtering Theory

2-1.
2-2.

2-3.

2-4.

The System.Model

Extrapolation of the Conditional Density Function
and Its Moments

Update of the Conditional Density Function and
Its Moments

The System Disturbance Model

Chapter III Discrete Filtering of a Linear Continuous System

3.1.
3-2.

3-3.

3-4.

Linear Filtering
State and Covariance Extrapolation

An Example of State Estimation of a Steam Turbo-
Generator System

The Problem of Filter Divergence

Chapter IV Approximate Filtering of Nonlinear Systems

4-1.
4-2.
4-3.
4-4.

NOnlinear Filtering
An Example of a NOnlinear Filtering Problem
Adaptive Filtering for NOnlinear Systems

A Discussion and an Example of Adaptive anlinear
Filtering

iii

10
11
15
16
18

24
28
33
33
37

44

49

Chapter V Conclusions
5-1. Thesis Review
S-2. Results and Conclusions

5-3. Suggestions for Further Research

LIST OF REFERENCES
APPENDIX A
APPENDIX B
APPENDIX C

iv

57

'57

58
60

61
69
72
74

Figure

Figure
Figure
Figure
Figure
Figure

Figure

3-1.
3-2.

4-1.

4-2.

4-3.

4-4.

4-5.

LIST OF FIGURES

The Block Diagram of a Steam Turbo-Generator System.

The RMS Estimation Error of Frequency Deviation
for the Steam.Turbo-Generator.

.Altitude Estimation Error for the First-Order
Filter of the Reentry Estimation Problem.

.Altitude Estimation Error for the Second-Order
Filter of the Reentry Estimation Problem.

.Altitude Estimation Error for the First-Order
Adaptive Filter of the Reentry Estimation Problem.

Altitude Estimation Error for the Second-Order
Adaptive Filter of the Reentry Estimation Problem.

RMS Altitude Estimation Error for Four Similar
Filters.

27

29

42

43

51

S4

I IAHRODUCTION

The subject of this thesis is the synthesis of discrete-time
filters for real-time estimation of the state of continuous-time
dynamic systems.

Historically, the theory of estimation dates back to the least-
squares orbit estimation of Gauss [GAUl] in 1809. More recently,
Wiener [WIEl] and Kblmogorov [KDLl] derived solutions for linear least-
squares estimation. 'Wiener studied continuous-time systems using
spectral methods to obtain minimum-variance unbiased estimates of the
scalar state by observing a signal where both the state and the obser-
vation are jointly wide-sense stationary and ergodic processes.
Kblmogorov used a recursive orthonormalization.procedure to solve the

analogous discrete-time problem.

In 1960, Kalman [KALI] published a recursive solution for the
minimum-variance unbiased estimate of the state of the Markov vector
process of a discrete-time system where the system is disturbed by a
random vector and the state is observed through a non-invertible trans-
formation. Kalman and Bucy [KAL2, KALS, BUCI] and Stratonovich [S'I'Rl]
extended this work to the continuous-time filter and incorporated the
problem of observations contaminated with additive random noise. In
addition to the filtering problem (estimating concurrently in time to

the observations), solutions were developed for the prediction problem

(estimating into the future) and the smoothing problem (estimating
after the time of the observations).

During the decade following the derivation of the Kalman-Bucy
filter, much effort was involved in further refinement and develOpment
of filtering theory. MDch of the motivation comes from optimal con-
trol theory where it was feund that a simple deterministic model would
no longer suffice for a realistic model of a system that is contami-
nated by noise and whose state can not be directly observed IBELI,
SARI]. The Kalman-Bucy filter was modified to serve in aero-space
applications for guidance, navigation, and orbit and reentry estimation
[SMIl, PURl]. In addition, applications in industrial process control,
electric power systems and system identification have been developed
[SAR1, CHAl]. More recently, applications in urban and freeway traffic
flow have received attention.

In a large class of optimal control problems, namely, for linear
systems with quadratic performance index and with Gaussian noise, the
optimal control law is a feedback law, which is a computed gain times
the state [BRY1, LEEl]. In this case, it can be shown that the com-
bined optimal estimation and control scheme is the optimal estimator
and the optimal controller found independently--this principle is known
as the separation principle [JOSl, POTl, WONl]. In the nonlinear case
and/or the non-quadratic performance index case, there is no known
separation theorem; nonetheless, the separation principle provides a
reasonable, approximate design approach.

Other approaches to the optimal control problem with inaccessible

states have been developed. In particular, the methods proposed by

Pearson [PEAl] and Ferguson [FERl] obviate the need for a filter. How-
ever, it should be noted that these particular methods are restricted
to optimal linear regulators and do not consider the case of noisy
measurements or system disturbances.

Linear filtering has been rather completely exploited. Since the
appearance of Kalman's original paper using the orthogonal properties
of random vectors to minimize the error criterion, other methods such
as least squares, maximum likelihood, minimum variance, calculus of
variations, and the maximum principle [BRY1, ATHl, H01] have been used
to prove the form of the optimal filter. Other problems such as non-
zero mean, cross-correlated, and auto-correlated system and measurement
noise have been considered [BRYZ]. Filter stability and convergence
have also been given considerable attention. In addition, the Kalman-
Bucy filter has been extended to the case of nonlinear systems, called
the extended Kalman filter, first-order filter, or the quasilinear
filter [FRIl, MOWl, SCHZ, H02].

The problems remaining in linear filtering are those that involve
mechanizing a filter that uses a suitable model, i.e., does not have
detrimental modeling errors, and that is feasible in terms of computa-
tional requirements.

Nonlinear filtering theory is not as completely developed as linear
filtering theory due to the inherently more complicated structure of
nonlinear systems where superposition does not hold. By using a
probabilistic approach, it is possible to derive the equations governing
the evolution of the conditional probability density function of the
state conditioned on the sequence of observations. It has been shown

that, in general, the optimal nonlinear filter requires an infinite

number of moments to realize the conditional probability density func-
tion [KUSl, KUSZ]. However, a finite approximation can be made that
may be satisfactory [KUSS, BUCZ].

~Another problem in nonlinear filtering is finding a global filter-
ing theory that will assure filter convergence. No general convergence
theory has been found that is equivalent to that in linear filtering
theory. Indeed, the problem is so difficult that proofs of filter
convergence are restricted to the scalar case without system distur-
bances [BIRl] . A

An approach is suggested in this thesis that improves the con-
vergence properties of both the linear and nonlinear filter whenever
modeling errors are prevalent. Previous work [DEN1, PINl, SMIl] has
been concerned with the problem of filter divergence. The idea of
adaptive filtering has been suggested by Jazwinski and Dennis [JAZZ,
DENZ], whereby an artificial system disturbance is constructed to
prevent filter divergence. In this thesis a new approach is suggested
to adapt the estimates to the measured statistics of the filtering
problem; namely, the filter gain matrix is corrected to provide con-
sistent filter statistics and measurement statistics.

In many problems, it is necessary to estimate the state of a
vector-valued continuous-time system. Often, it is not possible to
measure continuously, and even if continuous-time measurements are
available, only sampled data can be used to allow for the computation
time that is necessary even in high-speed digital computers. Thus, the
structure that emerges is a discrete-time filter for the continuous-

time system. This structure has some of the aspects of both continuous

problems and discrete problems.’ For the most part, the discrete-filter
problem is a sub-problem.of the discrete-filter continuous-system prob-
lem, whereas the continuous filter problem must account for the continu-
ous observation process.

One important problem involved in a continuous-time system is that
of modeling the random.processes. .Assuming the process to be continuous-
valued and composed of stationary independent increments leads naturally
to the Brownian motion process. It shall be assumed that this process
has a Gaussian distribution function with zero mean. In order to solve.
differential equations involving the Brownian motion process, the Ito
stochastic integral is employed.

To summarize the previous discussion, the purpose of this thesis
is to study the synthesis of filters. The discrete-time filter,
continuous-time dynamic system is emphasized. In Chapter 1, the general
filtering equations are developed, and some results are derived for simu-
lation of the stochastic processes. In Chapter 2, the linear filter
equations are given, some new computational methods are presented for
reducing computation'and increasing accuracy, these methods are applied
to a simulation example and finally, the problem of modeling errors is
discussed. In Chapter 4, a nonlinear filter is derived and used in a
simulation example. Finally, a new statistical method is developed that
adapts the nonlinear filter to modeling errors to provide consistent

elements for the nonlinear filter.

II. Analysis and Development of Filtering Theory

underlying the Subject of filtering theory is the theory of
probability, random variables, and stochastic processes. The central
principles of stochastic processes that are needed here are reviewed
in.Appendix.A. (See [FEL1, RAPl, D001], for example.) In this chapter,
the equations of filtering theory are developed for discrete-time

filtering of continuous-time, discrete-parameter, dynamic systems.

2-1. The System Model

 

A large class of dynamic systems can be modeled by the vector-
valued stochastic differential equation,
dx = f(x,t)dt + G(x,t)dB (2-1-1)

or its formal white-noise equivalent,

i = f(x,t) + G(x,t)W’ (2-1-2)
where x is an arrandom vector, 8 is an m-vector Brownian motion process,
f is a vector function, and G is a.matrix funCtion. The prOperties of
the Brownian motion process, its relation to the white-noise process
w(t), and its physical significance are given in [D001].

The existence and.uniqueness of the solution to (2-1-1) over an

interval T, namely,

x sz f(x,t)dt + fr G(x,t)dB (2-1-3)

(where p the integrals are interpreted in Appendix A) can be proven under

thelconditions that f and G are Lipschitz in x. (See [1101 or D001] .)

The Linear System Dynamic Model

 

In the case that f and G are not state dependent, the system model

is linear,

dx = A(t)x dt + B(t)dB, (2-1-4)
or

x = A(t)x + B(t)w (2-1-5)

where A is an n X n matrix, B is an n x m matrix, and B is a Brownian

motion process. The linear model is of great importance for several
reasons, one of which is that the solution can be written uniquely as

t
x(t) = acntpxcto) +ft accomodem (2-1-6)
0
where 6(t,to) is the transition matrix transforming the unforced
response over the interval [to,t] and the forced response is a stochas-
tic vector integral. If x(to) is taken as a Gaussian vector indepen-
dent of B, then x(t) is a Gaussian vector, and the covariance matrix

of the state P(t,s) can be written as

P(t,s) = E {x(t)x'(s)} (2-1-7)
= 6(t,to)P(to)d'(s,to) +
min(t,s) .
ft Manama-03' (or (mm-r
o
where min(t,s)
EiKS)B'(t)} =ft Q(r)dr . (2-1-8)

0

The Observation Model

 

The general form of the continuous observation model is the
stochastic differential equation

dy = h(x,t)dt + dv (2-1-9)
where y is an IL-dimensional observed quantity, h is a vector function

and v is a Brownian motion vector:
min(t,s)
E {v(t)v'(s)} = ft R('r)dt.

O

In the case of discrete-time measurements, the analogous model

of (2-1-9) is
y(th) = h(x(th), th) + v(th). ‘ (2-1-10)
The additive noise v is a white-noise sequence with covariance
B{v(th)v'(t£)} = R(th) 6 kt (2-1-11)

where 5M is the Kronecker delta function: 5a: l for ~12 = fl and
6M. = 0 for I2 f L In the case that h is not state dependent, (2-1-10)
may be written-as

y(th) = M(th)x(tk) + v(th) (2-1-12)
where “M is an IL x n matrix.

For the linear observation model and system model, a necessary
condition for estimating the state of the system is that the system
be observable [KALZ]. No equivalent theory exists for the general

nonlinear problem .

2-2. Extrapolation of the Conditional Density Function and Its Moments

 

The Markov, Brownian motion process x(t) satisfying the Ito

stochastic differential equation (2-1-1) has a probability density

function p(x,t) that is conditioned on the previous observations,
‘Y(tk) =‘{y(t1), ..., y(tk)}. Between measurements this conditional
probability density function, called the transitional probability

density function, p(x,t|Y(th)), is given by the Kolmogorov forward

 

 

equation,
, n n
épéiiil.= - g alp(x.t)fi]. §.z z azlpcx.t)[GQG']ij] 2-2-
t i=1 TXi i=1 i=1 ( 1)
BXiBXj

where fi is the ith_component of f and [GQG]ij is the ijth_element of

GQG'.. For a derivation and discussion of this expression and diffusion
processes in general, see [FELL 1101, DYNl, 1102, M310]. The initial
condition for (2-2-1) is the value at the prior measurement,

P(X:thl YCthD-

By the use of Ito‘s Theorem.A-3, the exact equations for the
extrapolation of the moments of p(x,t) between.measurements can be
determined. First, define the expected value of o in (A-ll) as

$CXCt)lth)) = E {¢(X)| Y(th)}=;(t(X) plx.t lY(tk)] dx (2-2-2)
Then the probability density, p(x,t|Y(th)) satisfies (2-2-1), and the

extrapolation of ¢ is

d$(x(t)lth) = wai-fmthn + §tr E {GQG'¢HIY(th)}dt,
th5t<th+1 (2-2-3)

from (A-ll) where each subscript denotes the partial derivative with
respect to x. The conditional moments are determined as follows:
for ¢(x) = x in (2-2-3), the conditional mean x(tlth), defined as

E{x(t)|th} is found from the solution of the differential equation,

10

d;(tlth) A

__..__.__. = E{f(X,t)lY(th)} = f(tlth):

dt
thf t < th+l' (2'2'4)

For ¢(x) = x x', the covariance matrix P(t|th), defined as
PCtlth) = E{[x - § (tlth)][x - Q (tlth)]'} (2-2-5)

is found from the solution of the matrix differential equation,

dP(tIth)

dt = E{Xf' ‘ X (tlth)f (tlth)}

+ E{fx' - f (tlth)x '(tlth)}

+ E{GQG'}, ths t <th+l (2-2-6)

Higher-order moments could be developed but will not be derived
here. The solution of (2-2-4) and (2-2-6) can only be approximated
for the general nonlinear estimation problem since an infinite number

of moments are needed to solve for p(x,t) [KUSl, KUSZ].

2-3. Update of the Conditional DensityFunction and Its Mbments

In the discrete-time filter, the other part of the filtering
problem is to compute the change in the conditional probability
density function.when additional information is available in the
form of a new measurement. This aspect of the problem is discussed
in various reports, [COXl, JAZl].

By the use of Bayes' rule, the probability density function

immediately after the measurement can be written as

p[x(th) Imp] = p[y(t,Q |x(t,,)1p[x(t,,) mtmn
pTY(tkllY(th_1)]

(2-3-1)

 

11

where

p[Y(th)|Y(th_1)] =](p[y(th)lx(th)1plxcthiIY(t,_1)1dx(th) (2-3-2)

and

P[Y(tk)lx(th)] P[Y(th) ’ h(x(th), th)]

= p[v(th)] (2-3~3)

which is a Gaussian distribution of zero mean and covariance given

by (2-1-11).

The expected value, x(tklthl is the conditional mean of the proba-
bility density of (2-3-1). MUltiplying by x and integrating over x
yields

;(tklth) = Eixct,)pty(t,)lxctk)1Itg} (2—3—4)

E{p[Y(tk)IX(tk)]|té}

 

where the superscript minus indicates a value just prior to the

measurement. The conditional covariance matrix P(th|th) is obtained

from the second moment of (2-3-1) and by expanding (2-2-5),

E{x(th)x'(th)p[y(tk)|x(th)]lté}
E{p[Y(tk3IX(th)]lté}

Pulzltlz) = (2-3-5)

 

The higher-order moments are not needed for this work and are
thus omitted. As was the case for the conditional moment extrapola-
tions, the updated moments depend on an approximation to solve for the
general first and second moments (2-3-4) and (2-3-5) since these
require the solution of the conditional probability density. The

approximate solutions are given in Section 4-1.

2-4. The System Disturbance Medel

 

In the case that the state model is of the form (2-1-5), then

the state covariance matrix can be written as

12

P(t) = 6(t,to)P(to)G'(t,to) +th’6(t,t)B(T)Q(t)B'(r)d'(t,r)dt
t

0
(2-4-1)
In this case, a discrete-time equation.may be used to model (2-1-5):

x(th+1) = d(tk+1,tk)x(th) + u(tk) (2-4-2)

This discrete representation can be justified by defining u(th) as a

white noise sequence with covariance matrix,

E{u(th) u'(tj)} =‘V(tk) 6h! (2-4-3)
where
tIt+1
V(th) = f C(th+1,T)B(T)Q(T)B' my (th+1,'r)d‘r (2-4-4)
t
h

A.numerical method of solving (2-4-4) is given in Section 3-2. For
purposes of digital computation, the simulation of (2-4-4) requires
a multivariate Gaussian distribution generator. The equations for
generating this distribution are given in.Appendix B.

.Although the above approach is straightforward, it avoids the
difficult paradoxes encountered by Bryson.and Ho [BRYl] and Bierman
[BIEl] in transforming from the continuous-time to the discrete-time
model.

Another method of approximating (2-1-5) is to write the
discrete-time equivalent as

x(t + At) & (I + A(t)At)x(t) +\/Zt B(t)u(th) (2-4-5)

where At is small and u(th) is a white-noise sequence with covariance

matrix Q(th), the same as the cov {w(t)}. The contribution of the

disturbance in (2-4-4) is

13

t + ‘
V(t ) =Ith 15(tk+1,T)B(T)Q(T)B'(1)5'(th+1,T)dT
h

i B(th)Q(tk)B'(th) At (2-4-6)
for tk+l - th+flt > 0.

The contribution of the disturbance in (2-4-5) is
V(th) = mm B(t)u(tk)u'(th)B'(t)\/At}
a B(tk)Q(tk)B'(th) At (2-4-7)

and therefore, for small enough value of At, (2-4-5) is equivalent at
each At interval to its continuous-time counterpart (2-1-5). This
approach may be considered where a sample function is desired with a
larger number of samples than in (2-4-2) and where the computation of
(2-4-4) is not necessary.

For a more general utilization of the property of discrete-time
modeling of continuous-time processes, consider (Z-l-l),

dx = f(x,t)dt + G(x,t)dB.
The discrete-time model is written as

x(t + At) = x(t) + f(x(t),t)At +\/A—t G(t)u(th) (2-4-8)

where u(th) is a white-noise sequence with covariance matrix Q(th)°
The contribution of the disturbance in (2-4-8) is the same as in
(2-4-7).

If the covariance matrix of the 8 process is given by (2-1-8),

then by Theorem.A-2,

E{ G(x,t)d3}= 0 (2-4-9)
At

14

and

E{j;tG(x,t)dB[ fAtG(x,t)dB]' =

I G(t,x)Q(t)G' (t,x)dt (2-4-10)
At

C(t.X)Q(t)G' (tut) At
for small At (i.e., the interval At must be small with respect to the
time rate of change of G and Q). Once again, the simulation generator

of this Gaussian multivariable disturbance is given in Appendix B.

III. Discrete Filtering of a Linear Continuous System

 

Since the conception of linear filtering theory in 1960, a great
amount of effort has been put forth, both in setting the theory straight
and in applying the theory to practice. 'Within a relatively short
period the theory was applied in space guidance and navigation programs
such as Ranger, Mariner, and Apollo. In addition, applications have
'been made of linear filtering theory to aircraft guidance, industrial
process cohtrol, and other detection problems. In most cases, it is
necessary to linearize the system differential equations about the
nominal state in order to apply a linear filter.

The explicit problems of deriving the optimal linear filter,
showing its stability, and proving asymptotic convergence are solved
[KAL1, KALZ, KALS, BUCl]. The problems remaining in linear filtering
theory are related to process modeling and filter mechanization.

The first part of this chapter is concerned.with developing the
equations for discrete-time filtering of the continuous-time system.
Next, some methods that have not appeared in the literature are
presented for mechanizing the linear filter. Then an example of an
application of linear filtering theory to a steam-turbine generator is
presented. Finally, a method for compensation for modeling errors is

presented.

15

16

3-1. Linear Filtering_

 

The linear set of differential equations that describe the linear
or linearized system may be written as

x = A(t)x + B(t)w I (3-1-1)
where A(t),is an integrable n x n matrix, B(t) is an integrable n X m
matrix and w is an m-vector white noise process with covariance
Q(t)6(t-r). .Although the white noise process is only a formal repre-
sentation of (2-1-4), its use is justified since in the solution of
(3-1-1), w(t)appears in the integral only and its integral, the
Brownian.motion process, is a well-defined function. The linear

?
observation process consists of sampled data of the form

V(tk) = MCtthCth) + V(tfz) . (3-1-2)
wherelM(th) is an.t X n linear non-invertifile transformation and v(th)
is an t-vector of a Gaussian white noise sequence with covariance
R(th). It is assumed that the system is observable.

The linear filtering problem is this: given initial estimates of
the state x(to) with Gaussian distribution and of P(to), the covariance
matrix, determine the mdnimum-variance recursive estimate of the system,
(3-1-1) and (3-1-2). Slightly different premises are assumed for the
discussion in Section 3-4.

With the modifications given in Section 2-4, the linear discrete-
time filter for the linear continuous-time system fits the framework
as derived by Kalman [KALl]. For the linear problem, the process x is
Gaussian and only the first and second moments of (2-2-1) and (2-3-1)
are needed, namely, the linear version of (2-2-4) and (2-2-6) for

estimation extrapolation and (2-3-4) and (2 3-5) for estimation update.

17

The linear filtering equations are given here without proof since
the derivation is commonly found in the literature and only a few steps
are required to reduce the above equations. Between measurements, the

state estimate and covariance matrix are

x(tkltk_1) a 0(th,th,1)x(th_1]th_1) (3-1-3)

P(tk|th_1) = 0(th,th_1)P(th_1Itk_1)¢'(tk,th_1)
+- fth w(th,t)'B(t)Q(t)B' (tw (thaw (3-1-4)
th_1
where o is the system transition.matrix. At the time of the measure-

ment, the state estimate and covariance matrix are
sent.) = ictklth_1) + Gummy.) - Mctkfiohltmn (3—1-5)
P(thlth) = P(th|tk-1) - G(th)M(th)P(thlth_1) (3-1-6)

where the filter gain is
ccth) = P(t,,lth_1)M'(tk)[M(tk)P(tk|th-1)M'(th) + Rctpl’l.
(3-1-7)
These five equations are the optimal linear filtering equations.
These equations are closely related to the maximum likelihood estimate
and the stochastic approximation method as shown by Ho [m1] .

When the system parameters are known exactly the solution of the
linear filtering equations is straightforward. However, the real-time
computer requirements may be exorbitant and unfeasible. The methods
discussed in the next section due to Frame and Needler [FRAZ] are
useful in decreasing computational time and increasing accuracy. In
the case that the system parameters are not known, the methods dis-

cussed in Section 3-4 can be used to correct for modeling errors.

18

3-2. State and Covariance Extrapolation

 

Minimizing computational time and maximizing computational
accuracy are crucial problems in.mechanizing an on-line filter for
filtering a complex process since these two criteria constrain
the filter capability in terms of (1) setting a lower limit on
sampling rate, (2) allowing local iterationsifor linearization of the
extended linear filter, and (3) allowing the system model to be suf-
ficiently complex to adequately model the physical system.

Two important aspects in both real-time system estimation and
control and in digital simulation of the same are extrapolation of the
state of the system.and extrapolation of the covariance matrix of the
system. In the case that the system may be described by (3-1-1) where
A and B are either time invariant or can be approximated as time
invariant over the time period of interest, then the computational
methods described in the sequel are useful in state and covariance
extrapolation between measurements.

The solution of (3-1-1) can be written as

t
x(t) = @(t-to)xo + Jr ¢(t-t)Bu(t)dr (3-2-1)
0
where the transition matrix ¢ is

<i>(t) = e“ = [235 Ahth/Iz! (3-2-2)

and to is taken to be zero.
Using the transition matrix to evaluate the transient response of
a linear time-invariant unforced system was described by Liou [L101]

and others [BRA1, MMSI]. Ganapathy and Rao [GAN1]* pointed out that

*For corrections to [GANl], see [RAOl].

19

the A.matrix has to be raised to successively higher powers to allow t
to be taken over a sufficiently large interval. They in turn show that
eAt can be expanded in an infinite series of scalars instead of matrices
and that.An'1_is the largest power of.A.that is required.

An alternative method requires slightly less computation and is
possibly more straightforward than [GANl]. This new method utilizes
the conjoint algorithm to compute the coefficients d1,...,dn of the

characteristic polynomial of A,

n-l +

n
A + dll . + d A + dn = 0 (3-2-3)

n-1

and the constant coefficient matrices, B0,...,Bn of the conjoint

matrix, defined by
n-1 n-2
B(A) = adjoint (A1 -.A) = B A + B A + ... B A+-B
° 1 n-2 n-1

(3-2-4)

where Bb = I and Bn = B_1 = 0. The conjoint algorithm described by

Frame [FRAl] yields the relationships,

Bk = AB}?1 + th (3-2-5)
db = -(1/h)tr.A Bh-l (3-2-6,
Define the transition matrix as
n-l
d>(t) = z B,z,(t) (3-2-7)
j=0 J J
where
0° j+h
z. t = z c 3-2-8
1( ) h ( )

 

12-0 (14-12)!

20

and the ck are determined so that

m
2: dc .=0form>0;c =d =1. (3-2-9
i=0 jM'j O O )

Then s1nce Bn=B_1=0, and zj==zj._1 forj>0wehave
, n
¢(t) - 2 B z (t) (3-2-10)

' i=0 1 1
n
= 2 AB z(t)+ )3 doIzj(t)
j=1 j-lj i=0 j
n d n 00 m
=AZB. z. t+I-—2 Zd-c .t/m! 3-2-11
ill 1"]. 1." ‘1( ) dt jao im=j J m'J ( )
or
5(t) =.A ¢(t) + I - o =.A ¢(t).l (3-2-12)

Also it should be noted that
n-l
<I>(0) = 2 BJ -jz (0) = BO 2 0(0) = B Gee; I (3-2-13)
j=0

Hence, we have the desired result:

¢(t) = eAt (3-2-14)

From equation (3-2-9) we solve for each Cl: using at most n multiplica-

tions:
m .
cm 121 djcm_, -j’ ms n
(3—2-15)
n
cm = - Z d c , m 2 n
j=1 Jm'J

By substituting these values of c into (3-2-8) and then solving

(2
(3-2-7), 9 (t) can be determined easily

21

The computational savings are considerable when compared with
using (3-2-2) truncated at some value b = N for N large compared to n.
On the other hand, when compared to the method utilizing the Hamilton-
Cayley theory [GANl], the savings in computation are less. In par-
ticular, the exponential series requires (N-l) n3 + N multiplications,

the Hamilton-Cayley method requires n4 - n3 - %n2 + %n + ZnN + N

multiplications, and the above method requires n4 - n3 - an + n

+ ZnN + N multiplications. The number of divisions is negligible and
additions are proportional to multiplications. Computer results con-
firmed the method and the previously mentioned computational savings;
for the example discussed later in this chapter, the computation was

about 30% as great as the computation of (3-2—2).

Covariance Extrapolation

 

In the case that u in (3-2-1) is a vector-valued, zero-mean,
white-noise process, then to determine the covariance matrix of the

error in x(t), it is necessary to evaluate*

t
V(t) =I ¢(t-t)Q<I>'(t-t)dr (3-2-16)
0

where

Q = E{B u(t)u' (t)B'}.

*This same quadratic form also appears in optimal control problems
with a quadratic performance index over a finite time interval. See

[1.1351] and [LEVI].

22

As pointed out by Levis [LEVI], errors are present from two sources:
truncation of the series expansion of ¢(t) in (3-2-2) and approximate
numerical integration.methods. In the following, the integral is
solved exactly and then a recursive relationship is derived to solve

for V'to the desired accuracy.

Theorem The solution of (3-2-16) where Q is a symmetric matrix is

 

given by
w t[2+1
V(t) = 2 G}, (3-2-17)
k=0 (h+l)!
where
Gh+1 = AG!2 + Gh'.A' (3-2-18)
and the value of Go is
G =Q. (3-2-19)

Proof: By substitution of variable, (3-2-16) can be rewritten as

V(t) = th @(t)Q¢'(r)dr (3-2-20)
0

After substituting (3-2-2) into (3-2-20) and integrating, the integral

 

becomes
m m . j+k+1
V(t) = )3 z AjQA'h t (3-2-21)
j=0 h=0 j!h!(j+h+l)
on k - _' [2+1
= z Z (kgAJQA',2 j t (3-2—22)

h=o i=0 1 (124-1)!

23

If G!2 is defined as

h .
(3,2 = z ( k)AjQA'M (3-2-23)
j=0
then
h1_.
GkA' = z (j)AjQA'M ,(3-2-24)
i=0

Replacing j by j-l in (3-2-23) it follows that
k+1h
Gk = 2 (j -1)A11QA'
j=l

k+1’1 (3-2-25)

and multiplying by.A gives

h+1 j h+1-j .
AG = 2 ('2 )A QA' (3-2-26)
k j= —l J 1

Thus, the sum of (3-2-24) and (3-2-26) is
Iz+l

G = (
h+1 j_0

h+1)AjQA'M 1'1 (3-2-27)

and by the symmetry property of Gk?

Gh+1 = AG,2 + (AGh)'. (3-2-28)

Thus, V(t) can be given exactly by

'0. tl2+1 -
V(t) = Z Gk
h=0 (h+1)!

 

(3-2-29)

where Gk is given by (3-2-25) and G0 = Q.
Although the benefits from utilizing this method cannot be evaluated
analytically, the increase in accuracy and the savings in computation as

compared to numerical integration are considerable in the case of large

24

t where it is necessary to have a large number of terms in the expansion,

as was realized in the following problemu *

3-3. An Example of State Estimation of a Steam Turbo-Generator System

 

A steam turbo-generator system is composed of a control system, a
steam boiler plant, and a turbo-generator. Disturbances occur to the
system in terms of a changing load caused by both changing consumer
power demands and changing power supplied by alternate sources. A
model can be constructed of the steam turbo-generator system by as-
suming that the coupling to other sources is relatively loose. See for
example, Kirchmayer [KIRl], Park [PAR1, PARZ], and Stanton [STAl]. Then
the system can be linearized around its operating point to obtain a
simulation model. The problem at hand is to estimate the state of the
system using this model. The model can be written in the form of
(3-1-1) and (3-1-2):

x = Ax + Bu

y = Mx + v
where in (3-1-1)

 

 

-1/T1 o o o Kl/Tl
1/T3 -1/T3 o o o
A = o /T3 -1/T o o ,
o KZ/i‘4 (1- eel/$14 -l/T4 o
L o -D/M l/M
0 ‘1
o
B = o ,
o
L 1 .1

 

 

*In the following example, less exact numerical integration required
approximately twice as much computation time as the above method.

25

and u is a random disturbance. The reference input is neglected for
purposes of this study. An explanation of the physical constants
appearing in the A matrix and M matrix as well as their numerical
values are given in Appendix C.

From a study of the random disturbances caused by changing loads,
it was determined that a random-walk.process would be an acceptable
model of this phenomenon. Since it can be shown [PAPI] that the
random-walk process is a limiting form of the Brownian motion process,
the Brownian motion process observed at discrete-time intervals will
exhibit the stochastic characteristics of the random-walk process.
Thus, in turn, the noise source can be pre-whitened by integrating a
white noise source to yield the Brownian motion process. The addition
of one state variable is necessary to account for the disturbance. The
state can be augmented with the differential equation,

16 =‘w

where W'is a white noise source. The augmented A matrix and B matrix

now become
”-1/T1 o o o icl/Tl T
1/T3 -l/T3 o o o

A = o l/T5 -l/T5 o o ,

o
o
o
K/T (1-K2)/T -l/T o o
2o 4 o4 -D/M4 l/M l/M
o o o o o

 

 

J

HOOOOO OCO

 

 

26

In Fig. 3-1, a block diagram.is shown of the overall system where each
component is represented by its Laplace Transform system function.

An acceptable model of the observation process depends on the
assumptions concerning the measurement devices. One choice is to
observe frequency deviation, x5, and frequency rate deviation, x5.
Thus, the observation matrix is
M=[0 o o l/M D/M 1m]

0 0 0 0 l 0
The noise source V(tk) can be modeled as a two-vector white-noise
sequence since the measurements are refined to the point that bias and
other gross errors are eliminated. A

The digital simulation'was performed using parameter data from
[PARl]. In addition, representative values were taken for the measure-
ment noise and system disturbance statistics. These values are given
in Appendix C.

The digital simulation of the steam generating plant including
the random effects, and the discrete-time filter demonstrated the
effectiveness of filter applications in problems of this type. By
using the filter as mechanized in Sections (3-1) and (3-2) it was pos-
sible to observe the estimation error time-history for various noise
sample functions. An improvement on simulating sample runs is to
study the estimation error dispersion by using the filter covariance
matrix of estimation error, starting with a large initial variance.

In this way, the need for performing Mente Carlo studies is obviated.

27

.ocwnnzu one nououocom onu we a

eofiuanw new 2.onuuoen mo peoEoe one new monsoooo nooHn Hocnw one .HHoEm ma osfio> Hocaeo: one

scum cofiuon>ou one an noocna mo wouoEnxouooo on xoe nonnon Emoum eeo_oono> eooum one .oenm>

Hmcneoe on» snow newuon>op xocoooohm one on ma .nofioonuuoo on .onon: Eoumxm one we moanofinm>

oumum.onu poo .mowconu omen genome an oomnoo 3 ooconnoumno eoecon onu .3 oocouomou women
eoHHoupcou one one moanowuo> one .Eoumxm nouonocou-onhoe eoopm m we Emumofim noon one .H-m .mnd

..L

 

 

 

 

n

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

$5135; ”5... .5 + _
m apex 1 . .x _ a
8.2425 . mamz: - moss. , 4 239m ,
235 mzaca. usammq 03mm .6528

 

.—|m

 

28

The tradeoff between improved filter performance and the following

factors can be evaluated: . I

l.~ The required sampling frequency which limits modeling complexity,
determines the maximum.computer cycle time, and determines measure-
ment transducer required sampling rates.

2. The quality and quantity of measurements that are required.

3. The effect of system disturbances.

Many simulation runs were performed while varying model parameters,
Such as sampling rate, the measurement accuracy and the magnitude of
the disturbance variance. A representative error-study time-history is
shown in Fig. 3-2.

In summary, the application of a Kalman filter to a steam-turbo
generator plant was simulated and the results were judged to be suc-
cessful. The methods proposed in Section (3-2) were used. The study
should be extended by evaluating the effect of (l) more refined model
development, (2) more exact noise and disturbance statistics, and
(3) application of a closed-loop optimal controller for minimizing
output frequency error. An alternative approach to improved process

modeling is presented in the next section.

3-4. The Problem of Filter Divergence

 

Early in the study of filtering applications, it was found that
although in theory the filter estimates converged asymptotically to the
state variables, in practice the filter estimates would diverge from

the state variables, especially where modeling errors were large relative

PER UNIT FREQUENCY ERROR x 10"

-6.0"‘

8.0- 7'9

 

7.0-

5.0'

 

4.0-

 

3.0‘

 

2.0-

 

 

' l l I
0 IO 20 30

TIME (SECONDS)

Fig. 3-2. The RMS Estimation Error of Frequency Deviation f or the
Steam Turbo-Generator. Starting with a very large error, the
RMS value converges relatively quickly to a steady-state error.

30

to system disturbance statistics. This outcome may be a consequence of
modeling error in any of the parameter values or covariance specifica-
tions [SCHl]. Since the filter has no measure of the consistency of
the estimation error covariance and the actual estimation error, it 15
possible for the covariance matrix to indicate a much more accurate
estimate than is actually being performed. The filter divergence is a
result of later observations being very lightly weighted since the
covariance matrix of the estimation error has converged to a smaller
quantity than warranted. .

Several approaches have been tried on the problem of filter
divergence. One approach is to augment the state vector with the model
parameters that are only approximately known [SMIl]. By estimating
these uncertain.parameters, this cause Of filter divergence is reduced.
The accompanying disadvantage is that a larger state vector leads to
greater computational requirements and_the possibility of establishing
an unobservable system. Another approach is to artificially include
a sufficiently large system disturbance covariance such that the filter
will not converge improperly [FITl]. Naturally, a good deal of modeling
and simulation needs to be performed to insert proper disturbance
statistics. One more approach suggested in [BERl] is to prevent the
determinant of the covariance matrix from decreasing below a given
reference value. A.parameter is chosen such that after the determinant
has converged for this length of time, the covariance is scaled so
that the determinant remains constant and the covariance matrix does
not go to zero.

With these methods, adjustments must be made strictly on the basis

ofTa-priori_knowledge obtained from measurements and simulation studies.

31

For this reason, methods are desired that are adaptive to the actual
estimation environment during the filtering process. 1

».As pointed out in [JAZZ], the filter divergence problem can be
judged by comparing the statistics of the actual measurement errors
with the covariance matrix of measurement errors. The only measure
of divergence available in the filtering problem is to compare the
actual measurement, y(th),'with the measurement estimate y(th) based
on the previous state estimate; namely, the measurement estimate
error e(tk) is A

B(th) = V(th) ‘ E{Y(th)|Y(th_1)} (3-4-1)

or 1

903k) = V(th) ' M'(tk)x(thlth-l)' ’ (3'4‘2)

The measurement estimate error covariance matrix 2(th'th-1) is defined

as

Z(th]th_l) = E{e(th)e'(th)} (3-4-3)

M(th)P(thIth_1)M'(th) + R(th) (3-4-4)

Since the mean measurement error estimate E{e(th)} is zero, it is pos-
sible to judge the consistency of the filter covariance with the actual
errors in the measurement estimate by comparing (3-4-4) with the statis-
tics of e(th)e'(th).

It is further shown in [JAZZ] that in the case of no system
disturbances and scalar measurements, a disturbance covariance can be
artificially created to adapt the identical noise inputs to the
statistics of the error in the measurement estimate. By maximizing

the probability density function p[e(th)] with respect to a scalar q

32

where qI represents the system disturbance, the most probable dis-
turbance is determined at each estimate. '

Another approach suggested by Andersen, et a1.,,[AND2] is to
obtain estimates of the state x, the transition matrix o, the input
covarianCe matrix V, and the noise covariance matrix R. The assumptions
are that the system is discrete, the parameters ¢, V, and R are constant,
the random sources are zero-mean, independent, and identically distri-
buted, but probably most importantly, the obéervation matrix M is‘
invertible, that is, each state variable is directly observable through
noise. .Although this latter restriction may be severe in many cases,
it should be noted that in a class of system applications this technique
may be valuable in taking adequate experimental measurements to ascer-
tain the true values of 0, V, and R and then using these values for
real-time estimation problems--assuming the other requirements are
satisfied.

A different approach is used in this thesis to handle the problem
of adapting the estimation procedure to the measurement data. This
proposed approach has applications in filtering linear systems with
parameter uncertainties; however, it has been used more extensively in
the non-linear filtering problem. Thus, the discussion of non-linear

adaptive filtering is deferred until the next chapter.

IV. Approximate Filtering of NOnlinear Systems

 

The subject of this chapter is discrete-time filtering of non—
linear continuous-time systems. The case’of system disturbances and
noisy measurements is covered. A second-order nonlinear filter is
developed, that is similar to the one suggested by Athans, et al.,
[ATHZ] but is extended to the case of system disturbances. Then a new
technique is developed by the author to adapt the filter to the measure-
ment statistics such that those statistics and the filter covariance

estimates are consistent.

4-1. Nonlinear Filtering_

 

The original linear filtering theory was adapted to nonlinear
systems by linearizing the system about the nominal value and using the
"first-order" or "quasilinear filter” as suggested in [FRIl, MOWI, SCHl]
during the early 60's.

By the middle 60's, the problem of nonlinear filtering was being
formulated [KUSl, KUSZ, BUCZ]. The exact equations for the propagation
of the state conditional probability density function conditioned on the
observation sequence are given by (2-2-1) and (2-3-1). However, an
infinite number of moments is required to realize the exact solution of
these equations. Thus, considerable effort has been extended to obtain

an approximate filtering solution to (2-2-1) and (2-3-1). These efforts

33

‘1!—

 

34

include [KUS3, SORl,.ATH2, JAZl].

In this chapter, the state model is assumed to be of the form,

1 = f(x,t) + G(x,t)w (4-1-1)
where the system dynamics function f is an n-vector which is twice-
differentiable with respect to x and once differentiable with respect
to t, the input matrix G(x,t) is a continuous matrix function of x,t
and as before, w is the zero-mean'white-noise disturbance with covari-
ance Q(t)6(t-T). The observation equation is of the discrete form,

V(th) = h(X(th)) + V(th) (4-1-2)

where the r vector h is also twice differentiable with respect to x
and once differentiable with respect to t, and as before, V(tk) is a
white-noise sequence with covariance R(th)5hg° And, finally, x(to) is
assumed to be a random vector with known Gaussian distribution.

Some definitions are needed to further expound the nonlinear filter.
The Jacobian matrix of f(x) is defined as A(x) where the ijth_entry of
A is
3f.

Ali-j. =3—Xj‘ ’ 4-3.1 = 19'.'°9n‘ (4-1-3)

Define the second partials of f as the Hessian matrix, F£(x) for
Z = l,...,n, where the ijth_entry of F3 is
32f

F£(x).. - —__£;_. 2,1 = 1,...,n. (4-1-4)

L1 - axiaxj

For the output function h(x,t), define the output Jacobian matrix

as M(x) where the ijth_entry is

af- .
o. = A ' = 000 = .00 4- -5
ML] 8x, 4' 1’ 3m, .1 1’ 9n ( 1 )

J

35

Define the second partials of h as the Hessian matrix, H£(x) for
2 = l,...,m, where the thh_entry of H1 is

32hz
H£(x),, = (4‘1‘6)
Lj Axiaxj

 

The function f may now be expanded in a Taylor series about the condi-
tional mean ;,

f(x) é f(x) + A(x)(x-x) + %.';1¢i(x-x)'Fi(x-;) ' (4-1-7)

&=
where ¢i is the basis vector, ¢t = [0...1...0]' with the one in the
ith_row. Likewise, the function h may be expanded to give

h(x) é h(i) + M(SE)(x-§) +%_— .314) p(x-QyHLu-i). (47-1-8)

g:
Similarly, G(x,t)Q(t)G'(x,t) can be expanded about x to give an approxi-
mation for the term E{GQG'}.

Terms that are higher than second order are neglected in these
expansions, although they could be included, at least in theory. In
practice, the error distributions are assumed to be Gaussian with zero
mean; thus all odd-numbered moments vanish [WOZl], and the fourth-order
moments may be decomposed into second-order moments [ANDl]. However, it
should be noted that the addition of higher-order terms does not greatly
affect the filtering equations for small estimation errors since they
converge to zero faster than the lower-order terms as the state estimate
converges to the state [SORl]. For larger errors, another approach

is suggested in Section 4-3.

36

For the nonlinear filter, the continuous-time dynamic system can
not be transformed to an equivalent discrete-time model as was the case
for the linear system by using (2-4-2) through (2-4-4). The differen-
tial equations for the extrapolation of the system state estimate must
be solved numerically.

The extrapolation equations are found by using the expansions for
f in (4-17) and for E{GQG'} in the moment equations (2-2-4) and (2-2-6):

1:: A In A (--
x f(x,t) + §T§1¢L tr[FL(x)P(t]tk_1], tk-ls t< t,2 4 l 9)

P(t|th_1) = A(x)P(tIth_1) + P(tlth_1)A'(x)
+ E{G(§)Q(t)c'(§)}, tk_1< t < th (4-1-10)

The updated equations for the nonlinear filter are found by using
the expansion for h in (4-1-8) in the moment equations (2-3-4) and
(2-3-5):

S‘ceklth) = icthltm) + G(tptyep - hole,» -

m A
'2’Lil‘bitrmfixfiknfltkIth)]] (4-1-11)
P(th|tk) = P(th|th_1) - G(th)M(x(th))P(th|th_1) (4-1-12)

where
ccth) = Pcthltmmv(ing)[Macthnpcthltflm'(Slap) I
+ R(tk) + 0(tk)]'1 (4-1-13)

37

The term 0(th) can be approximated as

_. 1 '" : ,
0(tk) - -4 221% tr[H£(x(th))P(tkItk_l)]
m ..
£51453 tr[H,-(X(th)‘)P(thlth_ID] (4-1-14)

when only first- and second-order terms are considered, or it may be

approximated as

Oij(th)_= %-tr [Hi(;(th))P(thlth_1)Hj(;(th))P(thlth_1)I (4’1‘15)

(0 is the £jth_entry of 0) when the fourth-order central moments are

41'
considered. These updated equations are derived in [ATH2, SORl, and
JAZl], only in slightly different form. The relative performance of
linear and nonlinear filters can only be judged by experimental results.
In the following section, an estimation problem is chosen that has
nonlinearities that are widely varied by making changes in parameter.

values. In this way, the value of the nonlinear filtering equations

are tested.

4-2. An Example of a anlinear Filtering Problem

 

In [WAGl, GRUl, BERl, WISl, ATHZ], the problem of reentry trajectory
estimation is discussed. Typically, both the reentry system dynamics
and the observation function are severely nonlinear and accurate esti-
mates are difficult to achieve. In this section this problem is dis-
cussed and the results are given for a particular example using first-

order and second-order filtering.

38

The Reentry Prdblem

 

It is assumed that the force of gravity is negligible when compared
to the aerodynamic effects. Thus, the system dynamics for vertical

motion of the reentry body can be written as [ALL1, ATHZ]

 

Sc =-x .

1 2 (4-2-1)
- _ CDAo 2

x2- 2m X2

where x1 is altitude, x2 is velocity, CD is the drag coefficient, A is
forward area, p is atmospheric density, and m is mass. The_atmospheric
density is assumed to be an exponential function of altitude,

o = po '2'“1 (4-2-2)

_ -1 .
where y = 5 x 10 5 ft . Since the ballistic parameter CD Ap/Zm 15 not
known it can be estimated by setting

X3 = C1) ADO/2m

and augmenting the state equations such that

X1 = 'XZ

x = -x2 x e’YXl (4-2-4)
2 3

x3 = 0.

It should be noted that no system disturbances are present in this
formulation.

The observation equation consists of a scalar radar range measure-
ment taken from an altitude H and a lateral displacement M. The discrete
measurements are contaminated with additive white noise so that the

expression for the observations is

39

.— 2 2 1/2 _
V(tk) " [M + (x1(t’2) ' H) ] + V(th) (4'2‘5)

where V(th) is assumed to be zero mean with covariance

' E{V(th)V(t£)} = ROM.

For simulation purposes, numerical values were given to the system
parameters. Although a great number of cases were tried, a typical
case with rather severe nonlinear effects was to offset the radar such
that the lateral range isiM = 100,000 feet and the radar altitude is
H - 100,000 feet. This evaluation produces the maximum measurement
nonlinearity simultaneously to the maximum deceleration nonlinearity.
The measurement uncertainty was chosen as R = 10,000 feetz. Measure-
ments were assumed to be taken at the rate of one per second.

The initial state vector was chosen as

x1(0) = 300,000 feet

x2(0) = 20,000 feet/second (4-2-6)

x3 (0) = 10'3 feet.1
Initial values for the state estimates were

§1(0) = 300,000 feet

A

x2(0)

20,000 feet/second (4-2-7)

5 1

3 x 10' feet-

523(0)

In addition, the initial estimates were assumed to be uncorrelated;

thus the initial state covariance matrix is diagonal with elements

p11(0) = 106 feet2
P22(0) = 4 x 106 feetz/second2 (4—2-8)
= 10‘4 feet'2

p33(0)

40

These particular values are identical to those used by Athans et a1.
[ATHZ] and.make a rather close comparison possible for computer simula-
tion studies. In addition, by varying the values of H and M as well as

the above initial values, a wide range of nonlinearities is effected.

Reentry Simulation Results

 

As stated in [ATHZ], "The validity and relative advantages and
disadvantages of --- (ad hoc assumptions for nonlinear filtering) ---
can only be evaluated and analyzed based on the evidence of real or
simulation results for particular systems." 'With this viewpoint, a
large number of simlations were performed for the nonlinear problem of
reentry estimation.

Since the estimation error depends on the particular random noise
sample function, it was necessary to perform a Mente Carlo study. As
opposed to the linear filtering case, the RMS error can not be assumed
to be identical to the square root of the appropriate covariance matrix

variance term since the nonlinear filter produces only approximate

 

minimum-variance estimates and covariance estimates. For this reason
several different statistics were found from the Monte Carlo study,
including the average of the estimation error, the average of the
absolute value of the estimation error, and the RMS estimation error,

defined as

N A
map = [11V ,zlcxcth) - x(tpfif” (4-2-9)
L:

where RMS(tk) is a vector. The value of N, the number of trials, was

chosen to be twenty.

41

The results of the Monte Carlo study using the previously given
initial conditions are depicted for the linear filter in Fig..4-l and
for the second-order filter in Fig. 4-2. For this particular simula-
tion, the absdlute value of the average error and the RMS error can be
Compared to results obtained by.Athans et al. [ATHZ].' Overall the I
results are similar, but it is interesting to note that whereas both
the average and.RMS final values of errors are within 10% of the
respective values obtained in [ATHZ] the average values have a much
smoother shape, indicating that more cancellation is taking place in
the averaging process than in [ATHB].

Several aspects of the Mente Carlo study of the linear and non-
linear filtering problem should be noted.

1. The second-order filter was more accurate than the first-order
filter, which tended to diverge as shown in the RMS curve of Fig.
4-1.

2. The effect of the second-order terms in the state estimate update
equation (4-l-ll) and the covariance update equation (4-1-12) were
nearly negligible whereas the state estimate extrapolation non-
linear terms (4-1-9) produced a significant improvement in the
estimation accuracy.

3. Both the linear and nonlinear filter accuracy appeared to be
virtually independent of the difference between computing with
regular or extended precision.

4. Although results are not shown for velocity x2, and ballistic para-

meter X3, the contrast between linear and nonlinear filtering
parallels that for altitude in that the second-order filter shows a

marked improvement in estimation accuracy.

(FEET)

IN POSITION ESTIMATE

ERROR

 

 

 

500- ‘
I I]: 42
a. _
II
: '1 .,
‘ . ,I II.
400- ”
I II
I If
I ',-\
I ,\=
I “
I ‘.
:3(>Clq I ‘:
I I“
I - r.
I I /\
I I /'\/' '-
, _ I‘ .- \
. . r’ ..
200- “ : \..—
I
I _ __ \_
‘I/\. I, \/ \\\
:5: ° . --
'00"l I, \\\/°. R
,1 \_I '\._\
I K...“
I N...—
I
0 AA I I 1
0 IO 20 30

TIME (SECONDS)

Fig. 4-1. Altitude Estimation Error for the First-Order Filter of
the Reentry Estimation Problem.
Absolute Value of Average Error
....... Average of the Absolute Value of Error
-.-.-.- Square Root of the Filter Altitude Variance

-..-..- RMS Error

(FEET)

IN POSITION ESTIMATE

ERROR

500 1

43

400 -1

300 -

ZOO-

 

lOO-I

 

 

 

 

 

O l I
0 IO 20 30

TIME (SECONDS)

Fig. 4-2. Altitude Estimation Error for the Second-Order Filter of
the Reentry Estimation Problem.

Absolute Value of Average Error

_______ Average of the Absolute Value of Error

....... Square Root of the Filter Altitude Variance
-.._.._ RMS Error

 

44

5. Neither filter exhibited an RMS error that converged as rapidly as
the corresponding square root of the diagonal term of the covariance
matrix, although the nonlinear filter reduced the discrepancy as

compared to the linear filter.

Other comparisons may be made but they are reserved until the next
sections where other factors of filtering theory and filter divergence

are discussed.

4-3. Adaptive Filtering for NOnlinear Systems

 

In the customary development of linear and nonlinear filtering, it
is assumed that the system parameters--state coefficients, initial state
estimation covariance, measurement coefficients, system disturbance
covariance, and measurement noise covariance--are known precisely for
the entire filtering time interval. In addition, it is assumed that
computation routines are exact, and finally, in the case of nonlinear
filtering, that the approximate filtering equations are accurate enough
that the Gaussian assumption for the estimation error distribution is
satisfactory.

In a large class of filtering applications these assumptions are
not warranted. The initial conditions may not be available or the
measurement noise and system disturbance covariances may change unex-
pectedly. And of course, the inaccuracies of numerical solutions and
using only approximate filtering must be considered.

It was shown in Section 3-4 that most approaches to the problem of
filter divergence are not on-line solutions--they depend on prior

information to the filtering process. Although the method suggested by

45

Jazwinski [JAZZ] is on-line, it applies to linear systems with scalar
measurements and without system disturbances. In addition, all of these
filter modifications have the effect of reducing the rate of convergence
since the state covariance matrix is a priori prevented from converging.

As reported in Denham and Pines [DENl] and Wishner [WISl], another
- approach to filtering a nonlinear system is to use the first-order
filter and perform iterations at the point of measurement, called the
iterated extended Kalman filter. Simulation results [WISl] have shown
that iterations are as effective as second-order filtering in reducing
divergence. Hewever, it should be pointed out that iterations increase
the computational load and limit the measurement sampling rate. For
this reason, in many systems, the better solution is to use the non-
iterative first- or second-order filter and to use a faster measurement
sampling rate. Thus, a requirement that is placed on divergence
correction techniques is that the computation load must not be signifi-
cantly increased.

As is the case in the linear filtering problem, Section 3-4, the
only measure of filter convergence performance is to compare the
statistics of the extrapolated measurement estimation error with the
covariance of the measurement estimation error. In the nonlinear case,
an approximation must be used to obtain the extrapolated measurement

estimation error,

ecth) yep - in.)
ecth) -- yep - M(§(t,,));<(tklth_1)- (4—3-1)

46

The measurement eStimate error covariance matrix 2(thltk'1) is defined

as

zcthlt ) -E{e(tk)e' (th ‘ ‘ (4-3-2)

k-l

= M(§(tk)P(thlthf1)M'(§(tk)) + R(th)

The statistics to be used for the extrapolated measurement estimation
error are defined to be of the form

N
8(th) = .20 aie(th-L)e'(th ) (4-3-3)

"4,

where N is a constant and “L is a weighting sequence on the outer
product of the past extrapolated measurement estimation errors such that

N
2 o- = 1.

i=0 i

In the case of ideal filtering, tr 5(th) should be comparable to
tr 2(tklth-l)’ that is, for Gaussian distribution of 5(th)’ tr 8(th)
should be less than tr 2(thltk-l) on an average of 68.2% of the
measurements since the diagonal elements of Z are variances. In the
case of filter divergence, 8(th) will be on the average large compared
to Z(th|tk_1) and the above will not hold.

New that this manifestation of filter divergence is realized, it
must be decided how this quantity is to be utilized to adapt the filter
to the real-time measurements.

The author's approach is to observe that the divergence of the

filter represents the fact that 2(thlth-l) has converged to a smaller

value than warranted by actual data, and this in turn means that the

47

state covariancelmatrix P(th|th_1) has converged to a smaller value than
warranted, which has the effect of diminishing the filter gain matrix
G(th). The discrepancy in the diverging filter can thus be alleviated
by using a gain correction.matrix Gc(tk) to compute the adaptive filter
gain Gash).

Ga(th) = Gc(th)G(tk) (4-3-4)

Then the updated adaptive state estimate is
xc(thlth) = x(thltk-l) + Ga(th)e(th) (4-3-5)

The gain correction matrix is selected such that the corrected

error

ec(th) = y(tk) - M xc(tk]th) (4-3-6)

is driven toward zero.

Consider the case of a scalar system. If the state variance has
converged incorrectly to a value P(th)/C(tk) where C(th) is the diver-
gence factor, then the gain G(th) is

act-,2) = P(t,,)M/(M2P(t,,) + cctpRcthn (4-3-7)
and the gain correction should be

M2P(th) + C(th)R(th)

 

Gc(tI2) = 2 (4-3-8)
M p t + R t
( I2) ( k)
A good measure of the divergence C is
S t
C(tlz) = ( (2) (4-3-9)

 

z t
( hltk'l)

48

:Thus, for the scalar system, the gain correction can be computed exactly
based on the measurement error statistic S(tk). However, it should be
'observed that the value of S(th) is dependent upon the measurement noise
sample function and that it is only an approximation to the true measure-
ment eétimation error variance. Therefore, a reasonable approximation

to (4-3-8) is

8(th)

N|I—I

(4-3-10)

 

Gc(th) =
2(thlt

h-l)

The general vector filtering problem.may be handled identically to
the above except the case of multivariable measurements requires sOme
modification. One modification is to reduce the measurements to scalar
measurements processed sequentially in time. Another approach that may
be considered is to use cross-correlation to identify those components
of the state vector that are correlated with the diverging measurements
estimates and apply the gain correction in accordance with the strength
of the cross-correlation.

An important property of the diverging filter is that the diver-
gence implies that the gain matrix is too small due to the unwarranted
convergence of the state matrix. Thus, by simply testing a criterion on
divergence and increasing the gain by an appropriate factor, such as
(4-3-10), the filter adapts to the measurement errors without the need

to increase the state covariance.

In order to assure filter stability, the gain correction is tested
to assure that oscillations do not occur. The test is performed by

examdning the corrected measurement estimation error,

49
ec(th) = y(tk) -:ch(th). (4-3-11)

If this error is of the same sign as e(th), then the adaptive
estimate is unaltered. However, if the corrected error is of the
opposite sign as e(tk), then the adaptive estimate has possibly been
over-corrected and to assure the filter will be stable, the adaptive
gain is reduced such that ec(th) = 0. Therefore, the adaptive filter

gain is

Ga(th) = Gc(th)G(th), e(tk)ec(tk) 3 0

]e(t )l ’ (4-3-12)

Ga(th) = h . Gc(th)G(th), e(th)ec(th) < o
16(tt11 * lec(th)l

Using this method of adapting the filter to the measurement esti-
mation error has the advantages that (1) it is not necessary to increase
the state covariance, which reduces or obviates filter convergence,

(2) it can be used even in the case that system disturbances are present,
and (3) the increased computational load is negligible with respect to

the normal amount of filtering computation.

1-4. .A Discussion and an Example of Adaptive N0nlinear Filtering

 

.A criterion to establish the performance of a filter is to compare
the RMS error of the estimates obtained by a.Monte Carlo study with
the square root of the variance obtained from the filter state covari-
ance matrix. But first, it is necessary to establish the validity of

the approximate filter covariance terms.

50

Although the accuracy of the covariance matrix cannot be explicitly
demonstrated, on the other hand it can be seen that the differential
equations governing its extrapolation are smoother than those of the
state estimate. This point has been supported by experimental evidence
in several ways. I
l. The covariance matrix does not vary significantly over the runs

used in the Mente Carlo study.
2. The covariance matrix does not vary regardless of whether the first-
or second-order filter is being used.
3. 'The covariance matrix is independent of the integration accuracy for
'widely varying integration intervals.
Any of these factors would greatly affect the state estimation error,
but the constancy of the covariance matrix supports the hypothesis that
the approximate filter covariance is nearly exact.
~ In the case of the second-order filter it can be seen from Figs.
4-1 and 4-2 for the reentry problem that the RMS error is much closer
to the covariance matrix values than is the case for the first-order
filter. NOnetheless the remaining discrepancy is still large, and the
purpose of using adaptive filtering in the reentry estimation problem
is to reduce this difference.

Fig. 4-3 shows the altitude estimation error for the reentry prob-
lem and is equivalent to Fig. 4-1 except the linear filter is replaced
by the linear adaptive filter. The results are significantly improved
over the linear filter; in fact, the RMS error is even less than the
RMS error for the second-order filter is shown in Fig. 4-2. Thus, in

this particular problem, the benefit of using this adaptive filter is

(FEET)

IN POSITION ESTIMATE

ERROR

500-

 

 

 

51
.I
A)
400 ' I '.
I
I)
‘ I
.‘-
300- I
I:
I
‘-
200- 1 A
I \.
I . \
\AS .. A.
7‘34 \‘ ,\\/‘~./\,
Ioo - 'CI’\ I \/ VAx/IR
I, \/ . \\ \/_
. x.‘. _
I I
0 IO 20
TIME (SECONDS)
Fig. 4-3. Altitude Estimation Error for the First-Order Adaptive

Filter of the Reentry Estimation Problem.

Absolute Value of Average Error

_______ Average of the Absolute Value of Error
_._._._ Square Root of the Filter Altitude Variance
-..,_.._ RMS Error

 

52

greater than that of using a second-order filter. Finally, the second-
order adaptive filter was mechanized and the equivalent altitude
errors are shown in Fig. 4-4. .As can be seen from the time-history of
the error curves, the discrepancy between.RMS and covariance-originated
terms for the adaptive second-order filter is reduced by a factor of
four as compared to the second-order filter in the latter part of the
run.

The four cases are summarized in the curves shown in Fig. 4-5. The
RMS altitude estimation error is shown for the linear, nonlinear,
adaptive, and non-adaptive filters as obtained from the Monte Carlo
studies. .

Several points should be noted with regard to these simulation
results.
1. .Although only x1, the altitude estimation error, is shown in Figs.

4-1 through 4-5, similar results were found for x2, the velocity

estimation error, and x the ballistic parameter, with regard to

3.
relative filter performance.

2. The value of S(th) for these particular results was based on a

A single sample of measurement error (i.e., N = l in (4-3-3)),
although similar results have been obtained for other trials, such
as an exponential weighting function oi in (4-3-3). The difference
in results was slight.

3. The effect of using extended precision and of including a greater

number of runs in the Monte Carlo study was negligible.

(FEET)

IN POSITION ESTIMATE

ERROR

5001

400 "I

300 -

200'

 

1‘; ¥
. "(0 -______.o

 

. " ‘l "
Ioo- It; ’ \\\’/\'

 

 

 

 

-‘I v \. \x
/ \7‘?‘ ..
I73}.-- '2
I 1 ”’41
0 IO 20 30
TIME (SECONDS)
Fig. 4-4. Altitude Estimation Error for the Second-Order Adaptive

Filter of the Reentry Estimation Problem.

 

Absolute Value of Average Error

-- Average of the Absolute Value of Error
._ Square Root of the Filter Altitude Variance
. ._ RMS Error

ERROR IN POSITION ESTIMATE (FEET)

5001

400‘

300‘

2001

l()C)‘

.OW] \?'<:;¢x\

 

Fig. 4-5.

~ 54

I
I y
- I
: I.
i F:
I: ’\
.- ' ’A \
. v .
. I o’\ -\\
‘\
\
A \
. I, , ‘\._..._._..
/\\ ..

 

I 1
IC) 2() 3C)

TIME (SECONDS)

RMS Altitude Estimation Error for Four Similar Filters.
First-order NOn-adaptive Filter

---- Second-order NOn-adaptive Filter

.-.- First-order Adaptive Filter
...... Second-order Adaptive Filter

55

4. The simulation results did not depend on the noise generator or on
the word length; the results obtained on the Michigan State Univer-
sity IBM 1800 computer were substantiated with further runs on the
Purdue University CDC 6500 computer.

In the above study, it should be noted that the assumed values of
initial state and initial state variance were consistent; that is, the
initial estimate of each state component was within one standard
deviation of the true state, also, the sample noise function was con-
sistent with the assumed noise variance, and finally, the integration
routine was virtually exact. Thus the only source for causing diver-
gence was the approximate filtering equations since the only inexactly
known parameter x3 was included in the augmented state. The other
possible uses of adaptive filtering become apparent when other system
inconsistencies are considered. For example, it was found that a much
coarser value of the integration interval could be used for the
adaptive filter without eliciting the problem of divergence that was
evident even in the second-order non-adaptive filter. Another example
is that even for the second-order filter it was found that the initial
state error could only deviate as much as three standard deviations of
the initial dispersion estimate for convergence to be assured and this
was for a relatively small initial state covariance matrix (4-2-8)
whereas for the adaptive filter a large initial error and covariance

matrix could be assumed without filter divergence.
Several aspects of adaptive filtering are worthy of further
examination. Using the statistic S(th) (4-3-3) leads to a very diffi-

cult question about what is the optimum weighting sequence

56

“2’ L = l,...,N and what is the optimum value of N. It appears that N
should be small compared to the number of measurements that are required
for the system to transcend a sizeable nonlinear range; in addition, for
heavily weighted past data, the effect of previous gain correction on
the state estimate should be considered. The proper function for the
adaptive filter gain (4-3-12) is not analytically determined and ap-
parently must be determined empirically.

In summary, the relative performance of the nonlinear filter and
the adaptive filter as they are developed in this chapter is shown in
Fig. 4-5 for state estimation of the atmospheric reentry problem and

the increase in accuracy for either filter is demonstrated while the

best performance is obtained by both nonlinear and adaptive filtering.

V. Conclusions.

 

The purpose of this thesis is the investigation of problems

involved in state estimation of a class of dynamic systems.

5-1. Thesis Review

 

The initial problems include determining a suitable mathematical
model for the system and stochastic model for the external disturbances.~
Due to the universality of the Brownian motion proces$ for modeling
real-world phenomena, as well as the availability of useful analytic
results, the stochastic disturbance model is assumed to be a Brownian
motion process (Appendix A and Section 2-1). Due to the large class of
dynamic systems that can be modeled by a set of ordinary differential
equations, the differential model (Z-l-l) is chosen for the dynamic
system model. The output of the system is assumed to be sampled
periodically as is the case in digital computer monitors and optimal
controllers; therefore, the discrete-time observation model (2-1-10)
is used. And finally, it is assumed that real-time estimates are
necessary as in the case of an on-line controller; thus the system
state estimator is assumed to be recursive. The results and conclu-
sions in this thesis are based on these assumptions.

Using the above assumptions, the filtering equations are develcped,
based on Kolmogorov's forward equation (2-2-1), Ito's Theorem A-3 aid
Bayes' rule (2-3-1). The worth of the filter mechanization is determined

57

58

by digital computer simulation; for this reason, the discrete-time
modeling of continuous-time stochastic processes (Section 2-4) is
derived.

In the case of linear filtering problems, two methods are pro-
posed for reducing computation time and increasing accuracy. An
example of state estimation of a linearized model of a steam turbo-
generator is given. And finally, methods are discussed for handling
the problem of filter divergence due to inexactly known models.

In the case of nonlinear filtering, the nonlinear filter is
developed and used to estimate the state of a reentry body. The

adaptive filter is suggested as a method to reduce the error in the

state estimate.

5-2. Results and Conclusions

The major results and conclusions of this thesis are based on over
200 simulations of the steam turbo-generator and the reentry body as well

as the results that have been referred to in the literature.*

1. The computational methods proposed in Section 3-2 for state and
state covariance matrix extrapolation can significantly reduce
computation time and increase accuracy in real-time State estimation.
The actual improvement depends on the size of the state and is more
significant with a larger dimensional system with more than two or

three state variables.

*In addition, some preliminary results on adaptive filtering of the Van
der Pol oscillator show favorable convergence properties.

59

The results of the digital computer simulation (Fig. 3-2) suggest
that real-time state estimation of a steam.turbo-generator is
feasible either for Optimal control of a single machine or for
area-wide control.

The second-order filter developed in Section 4-l successfully
estimated the state of a reentry body as determined by the digital
computer simulation. In all cases where the nonlinearities are
significant, the second-order filter has significantly less error
than the first-order filter. Examples of this are shown in Figs.
4-1 and 4-2 for Monte Carlo studies of reentry estimation.

The nonlinear filteris superior to the linear filter in the reentry
estimation simulation due to the second-order state estimate extra-
polation term; the second-order terms in~state update and‘covarianee
update have virtually no effect on filter performance.

The adaptive filter suggested in Section 4-3 successfully reduces
the estimation error significantly as compared to the second-order
filter for both linear and nonlinear adaptive filtering (Fig. 4-3
and 4-4). The fact that the adaptive second-order filter has an
RMS error that is comparable to the filter squarevrooted variance
indicates that within the given structure, any further improvements
in filter accuracy will be slight. These results are summarized in
Fig. 4-5.

The adaptive filter as developed in this thesis is an adaptive gain
filter; i.e., only the filter gain is adapted to the divergence

which is manifested in the extrapolated measurement estimation

60

error. Other methods of divergence compensation such as adaptive
noise covariance do not appear to be nearly as satisfactory as
adaptive gain filtering according to the simulation efforts of

this research effort.

5-3. Suggestions for Further Research

 

In the area of linear filtering, most of the research effort lies
in the area of applications and in filtering of inexactly'modeled
processes. Naturally, the problem of state estimation and concurrent
system identification are not resolved and require further definition.
(See Mbnahan.[MDNl] for a thorough discussion of system identification.)
An example is the problem of estimation and identification.where the
observation transformation is singular and the system description is
only partially known.

In the applications area, more data are required for further
refinement of the steam.turbo-generator state estimator presented in
Section 3-3. An obvious extension is to simulate cascading an optimal
controller with the optimal filter and to test the resulting improve-
ment in the quadratic performance index as compared to presently used
feedback regulation. Furthermore, the methods of Foshe and Elgerd
[F081] for inter-area Optimal control could be mechanized by coupling
the optimal controller with the optimal filter.

The second-order filter has been reviewed and its value has been
reviewed and its value has been demonstrated in this thesis. A.natural

question arises as to the utility of mechanizing even higher-order

61

filters. .Although this topic is certainly an area for future research,
the answer appears to be in the negative for several reasons. Higher-
order filtering requires higher—order derivatives, which in practice

may be difficult to compute either analytically or numerically. .Also,

higher-order filtering requires a great deal more computation (See

[SORl], for example), especially for a larger system. And finally,
higher-order terms appear to have a limited effect in filtering
convergence.

The use of adaptive filtering appears to be more promising for
improving filter performance of severely nonlinear systems, but this
conclusion should be tested over a wide range of examples and condi-
tions. In addition, the exact form of the optimal adaptive filter is

an open question and is the subject of further research.

LIST OF REFERENCES

[ALLl]

[ANDl]

[ANDZ]

[ATHl]

[ATHZ]

[BARl]

[BELl]

[BERl]

[3151]

[BIRl]

LIST OF REFERENCES

H.J. Allen and A.J. Eggers, Jr., ”A Study,of the Motion and
Aerodynamic Heating of Missiles Entering the Earth's Atmos-
phere at High Supersonic Speeds," NACA, washington, D.C.,
Tech. NOte 4047, October, 1957.

T.Wf Anderson, "An Introduction to Multivariate Statistical
Analysis," New York, Wiley, 1958._

W.N. Anderson, Jr., G.B. Kleindorfer, P.R. Kleindorfer and
‘M.B. Wbodroofe, "Consistent Estimates of the Parameters of
a Linear System," The Annals of Mathematical Statistics,
Vol. 40, No. 6 December, 1969, pp. 2064-2075.

‘M..Athans and E. Tse, UA.Direct Derivation of the Optimal
Linear Filter Using the Maximum Principle," IEEE Trans.
Automatic Control, Vol. AC-12, pp. 690-698, December, 1967.

 

 

M; Athans, R.P.'Wishner, and A. Bertolini, "Suboptimal State
EStimation for Continuous-time NOnlinear Systems from
Discrete NOisy Measurements, 1968 Joint Automatic Control
Conf., Ann.Arbor, Michigan, pp. 364-382, June, 1968.

M.S. Bartlett, Stochastic Processes, London, Cambridge
University Press, 1961.

 

R.B. Bellman, Adaptive Control Processes, pp. 129-131, 152-
158, Princeton University Press, Princeton, New Jersey, 1961.

 

.A. Bertolini, ”A Trajectory.Ana1ysis Program," Technical Note
1967-48, Lincoln Laboratory, MIT, Lexington, Mass., September,
1967.

G.J. Bierman, "A Remark Concerning Discrete Approximation to
White Neise and Covariance Calculation," IEEE Trans. Auto-
matic Control 14, 432-433, 1969.

 

 

‘M.W. Bird, “Asympotic Convergence of NOnlinear, Continuous-time
Filters," Ph.D. Thesis, Michigan State University, 1969.

62

[BRAII
[BRYl]

[BRYZ]
(BUC1]

[BUCZ]

[CHAl]

[COCl]

[cox1]
[DENl]

[DENZ]

[D001]
[DYNl]

[FELl]

[FERl]

63
F.H. Branin, Jr., "Computer Methods of Network Analysis,"
Proc. IEEE, VOl. 55, No. 11, pp. 1787-1801, NOvember, 1967.

A.E. Bryson, Jr. and Y.C. Ho, Optimization, Estimation and.
Control, Blaisdell Publishing Co., 1969.

 

.A.E. Bryson and D.E. Johansen, Linear Filtering for Time-
Varying Systems Using Measurements Containing Colored
NOise, IEEE Trans. Automatic Control 10, 4-10 (1965).

 

R.S. Bucy, "Optimum finite-time filters for a special non-
stationary class of inputs,” Johns Hopkins University, Appl.
Phys. Lab., Baltimore, Md., Internal Memo. BBB-600, 1959.

R.S. Bucy, "NOnlinear filtering theory, IEEE Trans. Automatic
Control, (Correspondence), vol. AC-10, p. 198, April, 1965.

 

S.Y. Chan and K. Chuang, ”A Study of Effects on Optimal Control
Systems Due to Variations in Plant Parameters,” Int. J.
Control, 1968, V01. 7, No. 3, pp. 281-298.

W.G. Cochran, and G.M. Cox, Experimental Designs, New York:
John Wiley and Sons, 1957.

 

H. Cox, On the Estimation of State Variables and Parameters
for Nbisy Dynamic Systems, IEEE Trans. Automatic Control 9,
5'12, 1964 o

 

W.F. Denham and S. Pines, Sequential Estimation when Measurement
Function Nonlinearity is Comparable to Measurement Error,
AIAA J. 4, 1071-1076 (1966).

A.R. Dennis, Functional Updating and Adaptive Noise Variance
Determination in Recursive-type Trajectory Estimators,
Special Projects Branch.Astrodynamics Conf. NASA-GSFC,
Greenbelt, Maryland, May, 1967.

 

J.L. Doob, Stochastic Processes. New York: Wiley, 1953.

 

E.B. Dynkin,.Markov Processes, Springer, Berlin, 1965.

 

W. Feller, An Introduction to Probability Theory and Its
Applications, Vbl. II, Wiley, New York, 1966.

 

 

J.D. Ferguson and Z.V. Rekasius, "Optimal Linear Control
Systems with Incomplete State Measurements," Proc. 6th Ann.
Allerton Conf. on Circuit and System Theory (Urbana, 111.,
1968), p. 135, April, 1969.

 

 

[F111]

[F081]

[FRAI]
[FRAZI

[FRIl]

[GANl]

[(34111]

[(39111]

[H31]

[H32]

[1101]
[11m]
[JAZl]

[JAZZ]

64

R. J. Fitzgerald, Error Divergence in Optimal Filtering Problems,‘
Second IFAC Syrup. Automatic Control Space, Vienna, Austria,
1967.

C. E. Fosha, 'Jr. and O. I. Elgerd, "The Megawatt- Frequency Control
Problan- -A New Approach Via (hatimal Control Theory," Dept. of
Elec. Engr., University of Florida, Gainesville, Florida.

J. S. Frame, "Matrix Functions and Applications," IEEE Spectrum,
Vol. 1, No. 6, pp. 123- 131, June, 1964.

J .3. Frame and M.A. Needler, "State and Covariance Matrix
Extrapolation," Proc. IE (to be published).

B. Friedland and I. Bernstein, "Estimation of the State of :a
nonlinear process in the presence of nongauss1an noise and
disturbances," J. Franklin Inst. , vol. 281, pp. 455-480, June,
1966.

 

S. Ganapathy and A. Suba Rao, "Transient Response Evaluation
from the State Transition Matrix," Proc. ‘IEEE, Vol. 57, No.
3, pp. 347-349, March, 1969.

F. K. Gauss, Theory of the Motion of the Heavenly Bodies Moving
about the Sin in Game Sections. New York: Dover Publica-
tions, Inc. , 1963 (reprint).

M. Gruber, "An Approach. to Target Tracking," M.I.T. Lincoln Lab. ,
Lexington, Mass. Tech. Note 1967-8, February, 1967.

 

Y. C. Ho, "On the Stochastic Approximation Method and Optimal
Filtering Theory," Journal of Mathematical Anaysis and
Applications 6, pp. 152- 154, 1962.

 

 

Y.C. Ho and R.C.K. Lee, "A Bayesian approach to problems in
stochastic estimation and control," IEEE Trans. Automatic
Control, vol. AC-9, pp. 333-339, October, 1964.

 

K. Ito, "Lectures on Stochastic Processes," Tata Institute for
Fundamental Research, Bombay, 1961.

K. Ito and H.P. McKean, Jr., Diffusion Processes and Their
Sample Paths, Academic Press, New York, 1964.

 

 

A.H. Jazwinski, Stochastic Processes and Filtering Theory,
Academic Press, New York, 1970.

 

A.H. Jazwinski, Adaptive Filtering, Proc. IFAC Symp. Multi-
variable Control Systems, Dusseldorf, Germany, OctOber, 1968,
2, 1-15.

 

 

[JOSl]

[KAIl]

[W1]

[KAI-2]

[W3]

[KIRl]

[KDLl]

[KUSl]

[KUSZ]

[KUS3]

[LEEl]

[LEEZ]

[LEVl]

R.D. Joseph and J.T. Tou, "On Linear Control Theory,

T.

R.

R.

.J . Kushner ,

.J. Kushner, "Nonlinear filtering:

6S

" AIEE
Trans. Applications and Industry, 80 (1961), pp. 193-196.

Kailath and P. Frost, 'Mathematical Modeling of Stochastic
Processes," Stochastic Problems in Control, 1968 Joint Auto-
matic Control Conference, Un1versity oFMichigan, Ann Arbor,
Mich., ASHE, pp. 1- 38, June, 1968.

 

E. Kalinan, "A new approach to linear filtering and prediction
problems," Trans. AM, J. Basic Engrg. , vol. 82, pp. 34- 45,
March, 1960.

 

E. Kalman, "New methods in Wiener filtering theory," Proc. >

lst %. on Engrg. Applications of Random Function Theory
0 a 1 1ty, I... Bogdanoff and F. Kozin, Eds. New York:

Wiley, 1963. , .

 

 

.E. Kalman and R.S. 'Bucy, "New results in linear filtering and

prediction theory," Trans. ASHE, J. Basic Engrg. , set. D vol.
83, pp. 95-107, Decemher, 1961.

 

.K. Kirchmayer, Economic Control of Interconnected Systems,

 

Wiley, New York, 1959.

.N. Kolmogorov, "Interpolation and extrapolation of stationary

random sequences," Bull. Acad. Sci. USSR, Math. Ser. vol. 5,
1941. A translationjias been puElished by the RAND Corp. ,
Santa Monica, Calif ., as Memo. RM-3090-PR.

 

.J. Kushner, "On the differential equations satisfied by

conditional probability densities of Markov processes, with
applications," J. SIAM Control, vol. 2, pp. 106-119, 1964.

 

"Dynamical equations for optimum nonlinear
filtering," J. Differential Equations, vol. 3, pp. 179-190,
1967.

 

the exact dynamical equa-
tions satisfied by the conditional mode," IEEE Trans. Auto-
matic Control, vol. AC-12, pp. 262-267, June, 1967.

 

 

.B. Lee and L. Markus, Foundations of Optimal Control Theory,

 

Wiley, New York, 1967.

.C.K. Lee, Optimal Estimation, Identification, and Control,

 

Research Monografii No. 28, M.I.T. Press, 1964.

.H. Levis, "Some Computational Aspects of the Matrix Exponen-

tial," IEEE Trans. of Automatic Control, Vol. AC-l4, No. 4,
pp. 410-411, August, 1969.

 

[LIOl]

[MASl]

[MCKl]

[MONl]

[MOWl]

[NAYl]

[PAPl]

[PARl]

[PARZ]

[PEAl]

[PINl]

[POTl]

[WM]

66

M.L. Liou, "A Nbvel Method of Evaluating Transient Response,"
Proc. IEEE, Vol. 54, No. 1, pp. 20-23, January, 1966.

E.J.'Mastascusa, "A Relation Between Lioufls Method and the
Fourth-order Runge-Kutta.Method for Evaluating Transient
Response," Proc. IEEE, V01. 57, No. 5, pp. 803-804, May, 1969.

H.P. McKean, Jr., Stochastic Integrals, Academic Press, New
York, 1969.

 

D.D. Monahan, "System Identification by Stochastic Approxima-
tion," Ph.D. Thesis, Michigan State University, 1970.

v.0. Mowery, "Least squares recursive differential correction
estimation in nonlinear problems,P IEEE Trans. Automatic
Control, v01. AC-9, pp. 399-407, OctOber,*I965.

 

T.H. Naylor, J.L. Balintty, D.S. Burdick, K. Chu,

Computer
Simulation Techniques, Wiley, New York, 1966.

 

A. Papoulis, Probability, Random Variables, and Stochastic
Processes, McGraw-Hill, New York, 1965.

 

G.L. Park, J. Preminger, A.P. WUrmlinger, and W. Panyan, "A
Deterministic and Probabilistic Analysis of Response of Power
Systems to Load Changes," Prog. Rt. 2, Div. Eng. Res.,
Michigan State University, East Lansing, Michigan, October,
1967.

G.L. Park, et al, "A Deterministic and Probabilistic Analysis
of Response of Power Systems to Load Changes," Final Report,
Div. of Engr. Research, Michigan State University, East
Lansing, Michigan, October, 1969.

J.B. Pearson and C.Y. Ding, "Compensator Design for Multi-
variable Linear Systems," IEEE Trans. Automatic Control,
Vol. AC-11, p. 130, April, 1969.

 

S. Pines, H. WOlf, A. Bailie, and J. Mohan, "Modifications of
the Goddard.Minimum Variance Program for the Processing of
Real Data," Analytica1.Mechanics Associates, Inc., Westbury,
New York, Rept. Contract NAS 5-2535, October, 1964.

J. E. Potter, "A Guidance- -Navigation Separation Theorem," Rep.
Re- 11, Experimental Astronomy Laboratory, Mass. Institute of
Technology, Cambridge, 1964.

P. E. Purser, M. A. Faget, N. F. Smith, Eds., Manned Spacecraft:
Engineering Design and_gperation Fairchild Publications,
Inc. , New York, 1965, pp. 226—231.

 

 

[RAOl]

'[SARll

[5cm]

[scnz]

[su11]

[son1]

[STAl]

[STRl]

[STRZ]

[WAGl]

[WIEl]

[WISl]

67"-

‘V.P.C. Rao and S. Ganapathy, "Comments on Transient Evaluation
from the State Transition Matrix," Proc. IEEE, Vol. 58, No.
5, p. 814, May, 1970. A

‘VIV.S. Sarma and B.L. Deekshatula, "Optimal Control When Some of
the State variables Are th Measurable," Int. J. Control, 1968,
V01. 7, N0. 3, pp. 251-256. .

 

F.H. Schlee, C.J. Standish and N.F. Toda, Divergence in the
Kalman Filter, AIAA.J. , 5, 1114-1120, 1967.

F.C. Schweppe, "Algorithms for estimating a reentry body's
position, velocity, and ballistic coefficient in real time
or for post flight analysis," M.I.T. Lincoln Lab., Lexington,
IMass., Group Rept. 1964-64, DDC609524, December, 1964.

G.L. Smith, S.F. Schmidt and L.A. McGee, "Application of
Statistical Filter Theory to the Optimal Estimation of Position
and Velocity on Board a Circunlunar Vehicle," NASA TR R-135,
1962.

H.W. Sorenson and A.R. Stubberud, "Recursive Filtering for
Systems with Small but Non-negligible Non-linearities,"
Int. J. Control, Vol. 7, N0. 3, 271-280, 1968. .

K.N. Stanton, "Estimation of Turbo-Alternator Transfer Functions
using Normal Operating Data," Proc. IEE, Vol. 112, p. 173,
1965.

R.L. Stratonovich, "Application of the theory of Markov process
f0r optimum filtration of signals," Radiotekhn. 1. Electron.,
vol. 5, no. 11, pp. 1751-1763, 1960. (English transl. in
Radio Engrg. and Electronics, vol. 1, pp. 1-19.)

 

 

R.L. Stratonovich, "A New Form of Representing Stochastic
Integrals and Equations, J. SIAM Control 4, 362-371, 1966.

W.E. wagner, "Re-entry Filtering, Prediction and Smoothing,"
J. Spacecraft Rockets 3, 1321-1327, 1966.

N. Wiener, The Extrapolation, Interpolation, and Smoothing_of
Stationary Time Series with Engineering Applications. New
York: Wiley, 1949. Originally issued as a classified report
by M}I.T. Radiation Lab., Cambridge, Mass., February, 1942.

 

R.P. Wishner, J.A. Tabaczynski and Mt Athans, "On the Estimation
of the State of Noisy Nonlinear Multivariable Systems," Proc.
IFAC Symp. on.Multivariab1e Control Systems, Dusseldorf,
Germany, October, 1968.

[WONl]

[MONZ]

[won]

68

NM. Wonham, "On the Separation Theorem of Stochastic Control,"
SIAM Journal on Control, May, 1963, Vol. 6, N0. 2.

 

E. Wong, and M. Zakai, "On the Relation Between Ordinary and
Stochastic Differential Equations, Internat. J. Engrg. Sci.
3, 213-229, 1965.

 

J .M. W02encraft and Irwin M. Jacobs, Princi les of Communication
Engineering, Wiley, New York, pp. 205-206, 1965.

 

APPENDICES

APPENDIX A

Stochastic Processes and Filtering Theory

 

In random processes, a Guassian white noise process serves as a
model for purely random effects that appear to be pure noise with
mutually independent components. However, care must be taken since,
by definition of its correlation function, a white-noise process is
almost nowhere continuous, and thus, not Riemann integrable. For this

reason, some properties of the Ito stochastic integral are needed [1101].

Definition A-l. Let B(t,w) be a Brownian motion process defined on

 

the interval te[to,tN] and let f(t,w) be measureable with respect to
B(t), the smallest. 0 field induced by (t,m), te[to,tN]. Then the Ito

stochastic integral is defined as

f tN N-l
f(t,w)d8(t) = 1.i.m. Z f(t»,w)[8(t. ,w)-B(t.,w)] (A-l)
tO 1+ 0 i=0 L L+1 L
where
f(t,w) = f(ti,w), ti 3 t < ti+1’ t0 3 t < tN QA-Z)
and 1 = max {t4+1 - ti}.

The need for careful attention to solving stochastic integral and dif-
ferential equations is elaborated on in Appendix A of [BIRl]. See also
[KAIl].

69

70

The following two important properties of the Ito stochastic inte-

gral shall be stated without proofs (see [ITOl] or [DOOl]).

Theorem A-2. Let the random vector f(t,w) satisfy the previous

definition, then

E{Jf f(t,w)d8(t)} = 0 ' (A-3)
T
and
I = . 1
E{ f1 f(t,w)dB(t)[ f1 f(t,w)d8(t)] PB f1" E{f(t,w) f (t,w)}dt
(A'4)
where P = cov{B(t,w)}.

B

These relations are needed to deal with stochastic differential
equations of the form I

1 = f(x,t) + G(x,t)é (A-S)
which.must be interpreted as a formal representation of the vector Ito
stochastic differential equation

dx = f(x,t)dt + G(x,t)dB (A-6)
The solution is of the form

x = IT f(x,t)dt + IT G(x,t)dB (N7)

where the first integral can be interpreted as an ordinary Riemann
integral but the second integral has meaning only if properly inter-

preted since B(t) is of unbounded variation.

71

An alternative to the earlier defined Ito stochastic integral is

the Stratonovich integral [STRZ]:

N'l x(ti ) + x(t

)
G( ,t)d8=l. Z £+1 , _
f X “:13 4?]. G( 2 9 t4) [B(ti'i‘l) 80-1)]

 

08-3)

where 1.i.m. is the mean-square limit. It has been shown by W0ng and
Zakai [WONZ] that if (A-6) is defined in terms of the Stratonovich

integral, then to be interpreted in the sense of Ito, it will become

Gx ~ ’
= f(x, t)dt+ +: §-§;-El- G(x,t)dt + G(x,t)dB (A-9)
where
BG(x, ____g_) n m 36. (x(t)
Gx ,t 2 2 I2 -1, , A-10

ax.
J

and the components of dB are independent. The difference between (A-6)
and (A-9) shows that the interpretation of the stochastic integral of
the stochastic differential equation used to model the physical system
governs the solution properties. The Ito stochastic integral will be
used exclusively in this work since it is more general; however, by
using the results of W0ng and Zakai, transformation can be made from the
Ito to the Stratonovich differential equation.

The main result from Ito [ITOl] is the following theorem.

Theorem A-3. (Ito) Given LA-6), then for ¢(x,t) a scalar-valued
function of x and t,

d¢ = ¢dt + o; dx + i'tr GQG'oxxdt (A-ll)
where the subscripts denote the partials, and the necessary partials do

exist.

APPENDIX B

The Mnltivariable Normal Distribution

 

If a Gaussian-distributed random vector x has expected value

:1 = E{x} (13-1)
and covariance

‘V = cov {x} = E{(x-x)(x-x)'} (B-Z)
then a simple method for generating random variates with the probability

density_function

= exp[-%(x-i)'V'lcx-i)1

 

MK) 1 2 (B-S)
[2w] /
is given in the sequel. (See also [NAYl] or [BARl].)
Using the transformation,
x = Tu + x (3’4)

Where i is an n-vector given by (B-l), T is a lower trangular matrix,
and u is an n-vector with unit-variance, zero-mean, independent,
Gaussian-distributed random components, then the variance of x can be
expressed as

‘V = TT'. (B-S)
Due to the symmetric form of V and the fact that T is triangular, T

can be determined from the relationships

72

73

11 1'1 V11 ‘
'-1
t.. = v.. a Z t. t. t.. 1 < ' < L < n B—6
41 (41 [2:1 4'2 1h” 11 J ' ( )
4-1
2 1/ .
t = - .
ii (vii k: tth) l < L < n

Thus, the generation of the multivariable random vector x with covari-
ance V requires the solution of (B-4) where T is determined in (B-6).
The components of u are generated independently by a Gaussian random

variate generator as in [COCl], for example. -

.APPENDIX C

Parameters of the Steam Turbo-Generator Filterinngroblem

The parameters of the steam plant model are listed below.
governor time constant

steam valve motor time constant
steam system time constant

reheat steam system time constant
reciprocal of regulation coefficient
reheat coefficient

effective rotary inertia of machine
damping torque coefficient

valve control signal

valve position deviation

torque output deviation

frequency deviation

load deviation

The values of the above constants are chosen for a typical machine:

~a ~a ~a ~a
LN
II I! II II

94 25 K1 = 62.24
3-8 K2 = 0.324
113-1 M = 18,221.6
4524.0 D a 5.94

74

75

All of these quantities are on a per unit basis and must be transformed
for real values. 2

In addition to the above parameters, values for the stochastic
model are chosen. The frequency and frequency rate measurement errors
are statistically independent with zero mean and variance
r11 = 0.0001
r22 = 0.0004
The system disturbance is a scalar with zero mean and variance
q = 0.0008
The initial value for the state covariance matrix is chosen to be very

large with uncorrelated components as follows:

p11(0) = 1.0
p22(0) = 1.0
p33(0) = 1.0
p44(0) = 1.0
p55(0) = 1.0
p66(0) = 1.0

The initial state estimate is chosen to be zero, §(0) = 0, whereas the
initial true state is given the value x(0) = 0.1; that is, each component
is set at 0.1. The sampling rate was varied widely in many different
simulation runs, but a value of 2 seconds was commonly used since this
value allows a reasonable computation interval and yet is not too large

with respect to the system dynamics.

mm“

m

Vldl
mm3
UHO

Iii

I‘WIIHHIH

3