
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00012992/00001
Material Information
 Title:
 An analysis of a new nonlinear estimation technique the statedependent Ricatti equation method
 Creator:
 Ewing, Craig M
 Publication Date:
 1999
 Language:
 English
 Physical Description:
 xi, 125 leaves : ill. ; 29 cm.
Subjects
 Subjects / Keywords:
 Acceleration ( jstor )
Aircraft maneuvers ( jstor ) Covariance ( jstor ) Density distributions ( jstor ) Noise measurement ( jstor ) Parameterization ( jstor ) Position errors ( jstor ) Simulations ( jstor ) Standard deviation ( jstor ) State estimation ( jstor ) Aerodynamic load  Mathematical models ( lcsh ) Aerospace Engineering, Mechanics, and Engineering Science thesis, Ph.D ( lcsh ) Dissertations, Academic  Aerospace Engineering, Mechanics, and Engineering Science  UF ( lcsh ) Nonlinear control theory ( lcsh ) Eglin AFB ( local )
 Genre:
 bibliography ( marcgt )
nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph.D.)University of Florida, 1999.
 Bibliography:
 Includes bibliographical references (leaves 122124).
 General Note:
 Typescript.
 General Note:
 Vita.
 Statement of Responsibility:
 by Craig M. Ewing.
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
 Resource Identifier:
 030472549 ( ALEPH )
43460574 ( OCLC )

Downloads 
This item has the following downloads:

Full Text 
AN ANALYSIS OF A NEW NONLINEAR ESTIMATION TECHNIQUE:
THE STATEDEPENDENT RICATTI EQUATION METHOD
By
CRAIG M. EWING
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
ACKNOWLEDGMENTS
I would like to thank my advisor, Dr. Norman FitzCoy, for his efforts in
continually pushing me to develop a better dissertation. I agree it is now something I can
look upon with pride. I also thank Dr. James Cloutier, the man behind the idea. Without
his help and guidance this work could never have been written. I would especially like to
thank my wife, Darsi, for her neverending support through these years of work. She
never let me give up.
TABLE OF CONTENTS
ACKN O W LED G M EN TS............................... .................................... ....................ii
L IST O F T A B LE S .............................................................................. .................... v
L IST O F FIG U R E S .................................................................... ..................................vi
A B ST R A C T ............................................... ................................................................... x
CHAPTERS
1 INTRODUCTION....................... ............................................. 1
2 LITERA TU R E SU R V EY ............................................................ ............................4
3 SCOPE OF W O RK .................................... ....................................................... 8
4 EXOATMOSPHERICGUIDANCEPROBLEM ENGAGEMENT SCENARIO ..... 10
5 EXOATMOSPHERICGUIDANCEPROBLEM FILTER DEVELOPMENT.......... 13
5.1 EK F D erivation ................................................................. ........................ 14
5.2 Bootstrap Estimator Derivation............................................... 19
5.3 SD REF D erivation ............................................................ ........................ 22
6 EXOATMOSPHERICGUIDANCEPROBLEM SIMULATION RESULTS........... 49
6.1 B aseline Perform ance.................. ...................... ................................................. 49
6.2 PDF Distribution Comparison................................................... 54
6.3 Sensitivity A analysis .............................................................. ..................... 64
6.4 SDREF Closed Loop Performance............................................... 79
6.5 Summary of Exoatmospheric Guidance Problem.............................................. 87
7 PENDULUM PROBLEM FILTER DEVELOPMENT.............................. .............. 88
7.1 SD REF D erivation .......................................... .......................................... 90
7.2 EKF D evelopm ent.............................................................. ....................... 91
8 PENDULUM PROBLEM SIMULATION RESULTS................................ ....92
iii
9 MSDREF DEVELOPMENT ................................................................. ............ 103
9.1 M SD REF D erivation........................................................................................103
9.2 Stability Proof................ .................................. .............................. 104
10 MSDREF PROBLEM SIMULATION RESULTS............................................. 112
11 CRITICAL ANALYSIS OF RESULTS ........................................ .................... 119
REFERENCES .................................................................................................... 122
BIOGRAPHICAL SKETCH................................................................................... 125
LIST OF TABLES
Table page
1. Filter simulation parameters baseline conditions............................. ............ 50
2. EKF and SDREF parameters onepercent initialization error.................................. 65
3. EKF and SDREF parameters fivepercent initialization error ................................ 69
4. EKF and SDREF measurement noise parameters level 1 ....................................... 72
5. EKF and SDREF measurement noise parameters level 2..................................... 72
6. M aneuvering target parameters................................... .............................. 76
7. Headon simulation parameters.......................... ....................... 81
8. 90degree beam shot simulation parameters ......................... .......... ........... .. 81
9. Tailchase simulation parameters.......................................... 82
10. Miss distance performance for varying scenarios................... ........................ 83
11. Miss distance performance for stressing target maneuvers.................................... 84
12. Miss distance performance for varying measurement noise................................. 85
13. Miss distance performance for varying initialization error...................................... 86
LIST OF FIGURES
Figure page
1. G generic intercept vehicle.......................................................... ...................... 11
2. Guidance system block diagram.......................................... 12
3. EKF position error with 3a standard deviations ................................................... 52
4. SDREF position error with 3a standard deviations ................................................52
5. Bootstrap position error with 3o standard deviations .......................... ............. 52
6. Comparison of EKF and SDREF covariance 3o standard deviations
for Y position .................................................................................................... 53
7. Expanded view of Figure 6 .................................................................................. 53
8. Posterior probability density functions at time = 1, 5, and 9 seconds
for X position .................................................................................................... 55
9. Posterior probability density function at time = 1, 5, and 9 seconds
for Y position ...... .............. ............................................................ 56
10. Monte Carlo posterior probability density functions at time = 1, 5, and 9
seconds for X position................................................................ ....................... 58
11. Monte Carlo posterior probability density functions at time = 1, 5, and 9
seconds for Y position............................................................... .......................... 59
12. EKF estimated and true acceleration............................................ ....................... 61
13. SDREF estimated and true acceleration......................................... ..................... 61
14. Bootstrap estimated and true acceleration..................... ......................61
15. Monte Carlo posterior probability density functions at time = 2, 2.5, and 3
seconds for thrust cutoff X position...................... ........................ 62
16. Monte Carlo posterior probability density functions at time = 2, 2.5, and 3
seconds for thrust cutoff Y position............................................ .................... 63
17. EKF position errors with 3o standard deviations for 1percent
initialization errors .................................................... .. .... ....................... 67
18. SDREF position errors with 3o standard deviations for 1percent
initialization errors ................................................................ .... ....................... 67
19. EKF position errors with 3a standard deviations for 5percent
initialization errors ................. .... .. ....................... 68
20. SDREF position errors with 3o standard deviations for 5percent
initialization errors ..................................................... .. ... ....................... 68
21. SDREF position errors with 3a standard deviations after new initialization
method, percent error .................................................................................... 70
22. SDREF velocity errors with 3O standard deviations after new initialization
method, 1percent error .... .......... ..... ....................... 70
23. SDREF acceleration errors with 3o standard deviations after new initialization
method, 1percent error ........................ ....................... 70
24. SDREF position errors with 3o standard deviations after new initialization
method, 5percent error................ ............................................................. 71
25. SDREF velocity errors with 3o standard deviations after new initialization
method, 5percent error ...................... ...................... 71
26. SDREF acceleration errors with 3a standard deviations after new initialization
method, 5percent error .............................................................. ...................... 71
27. EKF position errors with 3o7 standard deviations 1percent error in X,
100 microradian error in Y and Z............................................. ....................... 74
28. SDREF position errors with 3a standard deviations 1percent error in X,
100 microradian error in Y and Z.............................................. ....................... 74
29. EKF position errors with 3o standard deviations 10percent error in X,
500 microradian error in Y and Z............................................... ...................... 75
30. SDREF position errors with 3cr standard deviations 10percent error in X,
500 microradian error in Y and Z.............................................. ....................... 75
31. EKF position errors with 3a standard deviations for 10 deg/sec rotating
target m maneuver ................................................... ............... ...................... 77
32. SDREF position errors with 3a standard deviations for 10 deg/sec rotating
target m maneuver .......................................................... ................................ 77
33. EKF position errors with 3a standard deviations for dogleg target maneuver...... 78
34. SDREF position errors with 3u standard deviations for dogleg target maneuver. 78
35. A aspect angle description ........................................................ ............................ 80
36. Pendulum setup ............................................. ............................................... 88
37. EKF and SDREF estimates vs. true theta with no initialization error..................... 93
38. EKF and SDREF estimates vs. true thetadot with no initialization error .............. 93
39. EKF and SDREF estimation error for theta with no initialization error ............... 94
40. EKF and SDREF estimation error for thetadot with no initialization error............. 94
41. Monte Carlo probability density function for theta, with no initialization error,
tim e = 0.5 sec ........................................ .............. ................................... 95
42. Monte Carlo probability density function for thetadot, with no initialization
error, tim e = 0.5 sec......................... ... ............................ ....................... 95
43. Monte Carlo probability density function for theta, with no initialization error,
tim e = 2.0 sec ........................................ .................................................... 96
44. Monte Carlo probability density function for thetadot, with no initialization
error, tim e = 2.0 sec..................................... ............................................... 96
45. Monte Carlo probability density function for theta, with no initialization error,
time = 3.0 sec ................................................97
46. Monte Carlo probability density function for thetadot, with no initialization
error, tim e = 3.0 sec............................................................. ...................... 97
47. EKF and SDREF theta estimation error with initialization error............................ 99
48. EKF and SDREF thetadot estimation error with initialization error ..................... 99
49. Monte Carlo probability density function for theta, with initialization error,
time = 0.5 sec ................................................ 100
50. Monte Carlo probability density function for thetadot, with initialization
error, tim e = 0.5 sec............................................. ............... 100
51. Monte Carlo probability density function for theta, with initialization error,
tim e = 2 .0 sec ....................................................... ............. ......................... 10 1
52. Monte Carlo probability density function for thetadot, with initialization error,
tim e = 2 .0 sec ............................................................. ................................... 10 1
53. Monte Carlo probability density function for theta, with initialization error,
tim e = 3 .0 sec ...................................... ......................................................... 102
54. Monte Carlo probability density function for thetadot, with initialization error,
tim e = 3.0 sec ...................................................... .............. ......................... 102
55. True vs. estim ated X 1 state ....................................................................... 117
56. True vs. estimated X2 state ............................................................. 117
57. True vs. estimated X1 state, nonzero initial conditions....................................... 118
58. True vs. estimated X2 state, nonzero initial conditions.......................................118
Abstract of Dissertation Presented to the Graduate School
of The University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
AN ANALYSIS OF A NEW NONLINEAR ESTIMATION TECHNIQUE:
THE STATEDEPENDENT RICATTI EQUATION METHOD
By
Craig M. Ewing
August 1999
Chairman: Norman FitzCoy
Major Department: Aerospace Engineering, Mechanics, and Engineering Science
Research into nonlinear estimation techniques for terminal homing missiles has
been conducted for many decades. The terminal state estimator, also called the guidance
filter, is responsible for providing accurate estimates of target motion for use in guiding
the missile to a collision course with the target. Some form of the extendedKalman filter
(EKF) has become the standard estimation technique employed in most modem weapon
guidance systems. EKF linearization of nonlinear dynamics and/or measurements can
cause problems of divergence when confronted by highly nonlinear conditions. The
objective of this dissertation is to analyze a new nonlinear estimation technique that is
based on the parameterization of the nonlinearities. This parameterization converts the
nonlinear estimation problem into the form of a steadystate continuous Kalman filtering
problem with statedependent coefficients.
This new technique, called the statedependent Ricatti equation filter (SDREF),
allows the nonlinearities of the system to be fully incorporated into the filter design,
before stochastic uncertainties are imposed, without the need for linearization.
The SDREF was investigated in three problems: an exoatmospheric, terminal
homing, ballisticmissile intercept problem; a highly nonlinear pendulum example; and
an algorithmic loss of observability problem.
The exoatmospheric guidance problem examined nonlinear measurements with
linear dynamics. To investigate the SDREF when used with a combination of nonlinear
dynamics and nonlinear measurements, a highly nonlinear, twostate pendulum problem
was also examined.
While these problems were useful in gaining insight into the performance
characteristics of the SDREF, no formal proof of stability could be determined for the
original formulation of the estimator. The original SDREF solved an algebraic SDRE
that arose from an infinitetime horizon formulation of the nonlinear filtering problem. A
modification to the SDREF formulation was developed that led to a differential SDREF,
and a proof for local asymptotic stability was achieved. The modification removed the
infinitetime horizon assumption and integrated the coupled statedependent state and
covariance equations. This new form of the estimator is called the modified SDREF
(MSDREF). A problem involving algorithmic loss of observability was then examined.
This problem shows a performance advantage when using parameterization versus
linearization as in the EKF technique.
CHAPTER 1
INTRODUCTION
The extended Kalman filter (EKF) has become an accepted standard for both
tactical and strategic interceptor guidance problems. While some problems lend
themselves to this type of estimation strategy, systems with highly nonlinear dynamic
environments or misunderstood uncertainties can cause divergence of the filters. Many
alternative nonlinear estimation techniques have been researched in an attempt to
overcome this problem. The linearization of the nonlinear dynamics and/or
measurements is one area of concern. The EKF linearizes about the current estimate of
the state. If the state estimate is fairly accurate, this method has been shown to be
effective. If, however, the estimates are poor, linearization can cause errors in the filter,
which can lead to degraded performance or possible divergence of the state estimate.
It is this linearization that motivates the following research. This dissertation
investigates a new nonlinear estimation algorithm, which originates from a recent
nonlinear regulator design technique. It uses parameterization to bring the nonlinear
system to a linearlike structure with statedependent coefficients (SDC). The filtering
equivalent to the regulator technique can be realized in a recursive nonlinear estimation
algorithm, which has the structure of the steadystate continuous Kalman filter. The
Kalman gain is found by solving the statedependent Ricatti equation at each filter update
time. The filtering technique, called the statedependent Ricatti equation filter (SDREF),
2
allows the incorporation of nonlinear dynamics and/or measurements in the filter design
without the need to linearize the equations.
The SDREF was analyzed in a several ways to determine the benefits and issues
associated with this type of structure. By using a sampleddata Bayesian filter, known as
the bootstrap filter, the "true" probability density function (pdf) of the SDREF states was
monitored. This allowed an investigation into the differences in pdf evolution for
linearizing versus parameterizing filters. Finding the true pdf can be difficult; however
recent work, by the author on the bootstrap estimator, proved valuable for this analysis.
The bootstrap technique uses thousands of samples of the pdf along with Bayes' rule, to
estimate the true pdf of any given state. With this filter and a large number of samples,
the "true" state pdf can be closely approximated. An EKF was evaluated in a similar
manner for comparison purposes.
The SDREF was exercised in three problems. A realistic and detailed evaluation
of the filter was accomplished by examining an exoatmospheric guidance problem. The
SDREF, bootstrap filter, and EKF were coded into a sixdegreeoffreedom (6DOF)
FORTRAN simulation for analysis in an exoatmospheric terminal homing engagement
using passive and active sensor information. The sensor measurements included azimuth,
elevation and range. The characteristics of the SDREF and EKF were evaluated in
several different engagement scenarios. The SDREF was developed in a Cartesian
coordinate frame, which causes the measurements to be nonlinear, while the dynamics
remain linear.
The second problem was that of a highly nonlinear, twodimensional pendulum
problem. Evaluation of the SDREF when used with nonlinear dynamics and nonlinear
3
measurements was examined. The SDREF and EKF were coded into the visual
simulation environment VisSim. The performance of the filters and their pdf's were
compared. This academic problem was well suited to exploring any differences that
nonlinear dynamics may have on the performance or pdf characteristics of the SDREF.
No formal proof of stability could be determined using the above formulations of
the SDREF, so an alternate filter formulation was sought. The original SDREF was
based on an infinitetime horizon formulation of the nonlinear filtering problem that
resulted in an algebraic SDRE. A MSDREF was developed that utilizes a differential
SDRE. Because the Ricatti equation is state dependent, it must be integrated in
conjunction with the state estimates. A proof of local asymptotic stability was able to be
determined. The new MSDREF was examined in a problem consisting of loss of
algorithmic observability; and the results indicated a significant performance advantage
over the EKF.
This dissertation fully develops the requirements for the SDREF and MSDREF. It
includes an investigation of controllability and observability, which are pointwise and in
the linear sense and pertain to a valid SDRE solution. An exhaustive investigation of the
pdf characteristics, noise sensitivity, initialerror robustness, and closedloop performance
of the SDREF using the EKF as a comparative source was performed. A theoretical basis
for the MSDREF, when used as an observer, was shown. The following research
provides theoretical and practical information demonstrating that the SDREF and
MSDREF are viable candidates for nonlinear estimation, and for some problems can out
perform the EKF.
CHAPTER 2
LITERATURE SURVEY
The literature abounds with research of various nonlinear estimation techniques,
most of it mathematical and not application oriented. An elegant presentation of the
history of airtoair target state estimation techniques by Cloutier et al.' outlines a wide
variety of techniques, which include the Kalman filter, EKF, secondorder Gaussian filter
(SOGF), recursive nonlinear filter, recursive maximumlikelihood filter, and the assumed
density filter. These filters have also been examined in the context of single and
multiplemodel adaptive filters. Adaptive filters may make use of certain properties of
the residual or covariance to continuously change filter parameters to gain better filter
performance. On the other hand, multiplemodel filters use banks of filters to model
many possible target dynamic situations. Various methods are then used to either select
the results of one filter or combine that of several.
With the advent of the Strategic Defense Initiative (SDI), a new set of problems
arose when estimating the motion of a boosted intercontinental ballistic missile (ICBM)
with the use of passiveonly measurements. While the EKF is still the most widely
accepted estimator for this type of engagement, it has been shown to have degraded
performance when stressed by highly maneuverable targets.24 This is the same problem
that tactical filters experience in maintaining track of highly maneuverable aircraft or
missiles. When tracking a highly nonlinear maneuver, state estimates can become poor
and the linearization assumption begins to break down. The specific problems associated
4
with guidance filters for the exoatmospheric intercept problem have been addressed in
several research programs noted in the following paragraph.
Schmidt2 investigated an extension of work by Daum in which a closedform
solution to the FokkerPlanck equations for propagating the pdf was developed. A target
maneuverdetection algorithm was added to this formulation to enhance performance.
The emphasis in this research was to more accurately model the target acceleration based
on the rocket equation, instead of a Markovprocess. The target model caused the state
equations to become nonlinear as well. The results of the work did show an improvement
over the EKF technique; however this was accomplished, in part, through the use of the
maneuverdetection algorithm. Gido34 developed a nonlinear estimation algorithm,
which numerically implements the theoretical work of Kolmogorov. The algorithm
places a grid over the pdf and uses the FokkerPlanck equations coupled with Bayes' rule
to provide the state estimates. While this method showed improved performance over the
EKF, it is computationally intensive. In the author's masters' thesis at the University of
Florida,5 another novel nonlinear estimation technique was investigated. It involved
using thousands of samples of the pdf, along with the state dynamics and Bayes' rule to
propagate and update the pdf .7 The technique places no restrictions on the type of
statistics used to model noise. It also makes no linearization assumptions. While
computationally intensive, it is directly applicable to a massively parallel processing
architecture since the same operation is done to each sample.
While the bootstrap method may not be seen in any missile guidance system
design in the near future, it does provide a bruteforce method for obtaining a reasonable
representation of the true state PDF for a given problem. It is in this context that the filter
6
is used for the current research. While many of the above filters were shown to have
performance superior to that of the EKF, realistic implementation may be difficult.
Other techniques that have been proposed to solve the bearingsonly guidance
problem include the modified gain pseudomeasurement filter (MGPMF)8 and modified
gain EKF (MGEKF).9 Both use pseudomeasurements to avoid the need for linearization
around a "nominal" trajectory, such as is found in the EKF. These methods are restricted
to a special class of nonlinear functions, which can be manipulated into a pseudo
measurement form.
The construction of implementable nonlinear filters for nonlinear stochastic
systems remains a challenge. Numerous filtering algorithms, based upon series
expansions to approximately realize the conditional mean, have been suggested (see the
bibliography of Bucy et al.10) Once again, these techniques represented a heavy
computational burden, even for loworder dynamical systems.
The SDREF may alleviate some of the problems discussed so far. The technique
is derived from a new nonlinear regulator design methodology known as the state
dependent Ricatti equation (SDRE) control method.1112 The objective is to take the
nonlinear system described by
x= f(x)+g(x)u (21)
and transform it to the statedependent coefficient form
x= A(x)x + B(x)u (22)
where f (x)= A(x)x B(x)= g(x) (23)
7
using direct parameterization. Once in this form the SDREF is obtained by considering
the statedependent coefficients to be fixed and applying the linear continuous Kalman
filter to the "frozen" problem. The result is an algebraic SDRE whose structure is the
same as the steadystate Ricatti equation. The emphasis of this research is to develop and
analyze a guidance filter based on this technique. Initial work in the development of the
SDREF'3 demonstrates the filter in a simple twodimensional pendulum problem. Other
applications include induction machine estimation,14 simultaneous state and parameter
estimation, 5 and gyroless satellite angular rate estimation.16 The current research
extends their work to include an exhaustive investigation into the internal properties of
this filter in a strategic interceptor guidance problem, as well as in the frictionpendulum
problem. Insight into the internal representation of the pdf when using parameterization
versus linearization is examined, as well as the necessary conditions for stability. The
first example presented here uses only nonlinear measurements, while the pendulum
problem discussed herein adds the complexity of nonlinear dynamics as well. A third
example demonstrates the unique parameterization property of the SDREF when
confronted with a loss of algorithmic observability.
CHAPTER 3
SCOPE OF WORK
This chapter outlines the scope of the work in this dissertation. Throughout the
design and analysis of the SDREF, the EKF was used to compare performance issues
between linearized and parameterized methods. Similarly a sampleddata bootstrap filter
was used to compare probability distribution functions between the different nonlinear
filters in the exoatmospheric guidance problem. The three filters are theoretically derived
and applicationspecific equations are formulated for each problem. Controllability and
observability requirements for the SDREF algorithm are analyzed and determined.
Several performance measures are investigated that include sensitivity to initial errors,
noise parameters, and target maneuvers. In addition, a somewhat nonstandard
methodology, the bootstrap filter, is used to investigate the generally unobservable
evolution of the pdf within the filter. The pdf was examined at distinct time steps to
determine whether the parameterized filter pdf is different from the linearized filter. All
three filters were coded in FORTRAN and installed in a 6DOF simulation for evaluation
in an exoatmospheric terminal homing engagement.
To examine the effects of both dynamic and measurement nonlinearities, the
SDREF and EKF were also coded in VisSim (a visual simulation tool for a simple, but
highly nonlinear, pendulum problem). Performance was measured in a manner similar to
the exoatmospheric guidance problem. Because no formal proof for stability could be
found for the SDREF, modifications were made that allowed a detailed proof of local
8
9
asymptotic stability to be found. The MSDREF is fully derived and an examination of its
performance against the EKF for a problem consisting of loss of algorithmic observability
is also explored.
CHAPTER 4
EXOATMOSPHERICGUIDANCEPROBLEM ENGAGEMENT SCENARIO
The filters were initially analyzed in the context of an exoatmospheric terminal
homing engagement of an interceptor tracking a thrusting ICBM for approximately the
last 10 seconds before impact. The interceptor was closing on the target at approximately
10 kilometers per second. The target may perform thrust cutoff or other evasive
maneuvers, which provide for a stressing nonlinear engagement scenario in which to test
the filters. At the start of the engagement, the guidance filter was given erroneous initial
state and covariance estimates. The sensor is assumed to have acquired the target and the
guidance filter is ready to estimate the relative motion between the two vehicles using
passiveonly measurements. The estimates were used by a standard augmented
proportionalnavigation (PRONAV) guidance law to command accelerations that keep
the interceptor on a collision course with the target. No errors associated with the attitude
control system were used as this would complicate the analysis and cloud the
performance comparison.
Figure 1 shows a model of the interceptor, a generic 4kilogram vehicle
configured with divert and attitude control thrusters. A generic, SS18like ICBM model
was used to generate the target motion. The interceptor may have many error sources,
some of which include handover errors, pointing errors, thrust errors, centerofgravity
location errors, and seeker noise. Initially many of these noise sources were not used so
that the filter characteristics could be more easily discerned. Several engagement closing
10
11
angles were examined to allow for more stressing conditions (i.e., tail chase, noseon,
beam shot). Once proper filter operation was assured, noise sources were added to more
realistically emulate the true intercept scenario, as well as to identify weaknesses in filter
performance.
Back View
SDivert Thrusters
Top View Z
Divert Thrusters Back View
Attitude Control Thrusters
X Z Z
Figure 1. Generic intercept vehicle.
The simulation used was the Guided Simulation Program Exoatmospheric (GSP
EXO).' It is a FORTRAN simulation developed by the General Electric Company for
the Air Force. Figure 2 shows a block diagram representation of the subsystems
modeled.
Figure 2. Guidance system block diagram
To allow for rapid and efficient evaluation of filtering concepts, the author has
modified the guidance system. It has a MonteCarlo capability with an outer loop around
the MonteCarlo loop, to allow for parameter trades. Analysis was performed using a
personal computer; however, analysis can be accomplished on other platforms, including
the SUN, VAX, and CRAY.
CHAPTER 5
EXOATMOSPHERICGUIDANCEPROBLEM FILTER DEVELOPMENT
The following paragraphs outline the formulations of the EKF, bootstrap filter,
SDREF, and MSDREF that were used for this research. A brief theoretical description of
each filter is given, followed by the applicationspecific equations.
There are three filter types under consideration in this research paper. The
mathematical derivations for all of the estimators are fully developed and all assumptions
regarding noise characteristics and linearization are explained. The filters are presented
in coordination with the problems that were used in examining their performance. The
SDREF theory and derivation used in the exoatmospheric problem is shown in this
Chapter, the SDREF pendulum problem and derivation is presented in Chapter 7, while
the MSDREF stability proof and analysis is found in Chapter 8.
Many variations of SDREF are possible due to the infinite number of
parameterizations available. Initially, simple parameterizations were used. Another
important factor in evaluating performance differences between the three filters was the
selection of the noise characteristics. The EKF performance may be enhanced in some
cases by artistically "tuning" the noise parameters within the filter. To assure an unbiased
comparison between the filters, noise parameters were set to the same nominal values for
all filters. Due to several limitations on the bootstrap filter, developed in a thesis by this
author, it was used only for pdf comparisons in the analysis of the strategic guidance
problem. The bootstrap filter needs further development to be useful in stressing
13
14
engagements as the number of samples required for adequate performance makes it
intractable to run for large numbers of comparison.
The SDREF and EKF were analyzed in great detail to assess key performance
characteristics including convergence, stability, and state estimation error. A key
performance critique was the comparison of the evolution of the SDREF and EKF pdf's
to that obtained by using the bootstrap filter. Interesting insight into the differences
between how the SDREF and EKF model nonlinearities was observed.
There are many issues involved with the SDREF development, including correct
parameterization of the nonlinear measurement equations, investigation into the
requirements for solution of the SDREF, and the stability of the SDREF. These areas
were explored to determine an accurate set of state and measurement equations and a
feasible method to solve them.
An initial non stressing scenario was selected to establish a benign comparative
baseline performance assessment. The pdfs were compared using this baseline scenario
as well as a nonlinear maneuvering target scenario. Once the baseline was established,
the filters were tested against more stressing engagements and intercept angles.
Performance was measured by examining the state estimate error and 3sigma standard
deviations, as well as comparing the state estimate to the true state. The author has found
that this is an accurate method of evaluating filter performance.
5.1 EKF Derivation
The derivation for the EKF can be found in many sources,1819 and is presented
here for completeness. The dynamics for the intercept problem at hand are linear, while
the measurement equations are nonlinear. For this reason a combination of linear and
EKF equations were used and implemented in discrete form. Given a linear, timevarying
differential equation for a dynamic system
x(t) = F(t)x(t) + G(t)u(t) + L(t)w(t) (5.1)
where x(0) is given, u(t) is a deterministic input, and w(t) is a random disturbance, the
corresponding discrete time solution is
t,
x(t) = D(t, t_,)x(t_) + fQ(tr .,)[G()u(r) + L(r)w(r)dr. (5.2)
1.
If u and w are assumed piecewise constant over the sampling interval, this can be
written as
Xk = ~ kXkI + Tkluk +Ak ,_,wk (5.3)
where rk_ k = J (tk,r)G(T)udr (5.4)
and AkIW, =j (tk,r)L(r)wdr. (5.5)
1ti
However if F, G, and L are also considered constant over the sampling interval, AT, then
4(At) = e" (5.6)
F(A) = D(At)[l, D'(At)]F'G (5.7)
A(Ar)= 'D(At)[,, '(At)]F'L. (5.8)
Assuming then that
E{xo) =m, E{(xxo)(xxo)T) = P
E{w)= 0 E{wwwkT)=Q; and
E{ww,_,} =0
the state propagation equation becomes
xk = (), xk + rk,Uk_ + A,_(0) (5.9)
The covariance matrix can be propagated in a similar manner by
P. = l,P k, l'r+ k (5.10)
with Qki = AkQAk. (5.11)
To relate the covariance matrix of random sequences to the spectral density of the random
process for the continuous disturbance input w(t) we note that
E{w(t)wT ()}= Q'c (t)S(t z) (5.12)
The equation for Qk,_then becomes
Qk1= i (tk')L(r)Q'c ( r)r(t t)dr (5.13)
k;
This now represents the integrated effect of the continuous disturbance input w(t).
Thus the two equations necessary for the propagation of the state and covariance
are
xk = kl,_l + k_Uk_, (5.14)
P'k = kIPk_ ITi' + Qki (5.15)
where Ok1 is defined by eFM and 1_, is define in equation (5.7).
The measurement process for the filter is now determined. Given the continuous,
nonlinear measurement equation
z(t) = h(x(t))+ v(t). (5.16)
The discrete version can be written as
z, = h(xk(t))+v (5.17)
17
A complete derivation of the EKF update equations can be found in,18 however the final
solution for the calculation of the Kalman gain, and the update of the covariance and state
are
k (+)= ()+ K [zk hk(k())] (5.18)
K, = 'k()HkT (i())[Hk(ik())Pk ()H'(k (())+Rk (5.19)
Pk(+) = [I KkH,(())]Pk()[ KkHk (k ())I + KRKkT (5.20)
where ik () and Pk ()come propagation equations (5.14) and (5.15) and
Hk (xk = (5.21)
with E{v } =0 and E {vkvk } = Rk. Now (5.14) and (5.15) and (5.18) through (5.20)
define the full set of EKF equations.
For the exoatmospheric intercept problem the state consists of
sxR
sYR
SZR
xR
xk = YR (5.22)
VZR
axr
ay,
az,
where the states are relative position, relative velocity, and target acceleration. The
system equations are
sR =VR (5.23)
vR = AT A,, (5.24)
ar, = Aa + w (5.25)
Putting this in state space notation,
Ss 1 0 S 0 0 0 000
S1=0 I0 v, + 00 I 0 a co 0 + 0 0 (5.26)
a, 0 0 AJla, 0 0 o0 0 0 L o wr
From this equation the state transition matrix is computed as
I IAt I(e "+ At 1)/A'
0(At)= 0 1 I(le)/IA (5.27)
0 0 le'
where I is the 3x3 identity matrix.
Assuming the commanded acceleration remains constant over the filtering interval
IAt2A ,I
2
SD(tk,,r)G(r)udT= AtAcoA (5.28)
0
The state propagation equation becomes
At AcaI
+ tAO1 (5.29)
0
The covariance propagation proceeds in a similar fashion. The covariance can then be
propagated by
P = [ <]_, []T + Q. (5.30)
The measurement equation is
zk = hk (,(tk) + v,
(5.31)
Range
where h (i(tk))= zE / Range Range = x + + R (5.32)
y, I Range
and vk = N(O, R) (5.33)
is the Gaussian (normal), random vector with zero mean, and covariance Rk.
This completes the derivation and equation formulation for the ninestate EKF
needed for the interceptor guidance system. Given an initial estimate of the state x(0)
and the covariance matrix P(0) equations (5.14) and (5.15) and (5.18) through (5.20) can
be solved in a recursive fashion to achieve an estimate of the relative vehicle motion.
5.2 Bootstrap Estimator Derivation
The bootstrap filter is derived using the recursive Bayesian estimation
relationships. To begin the discrete dynamic state equation is given by
Xk =fk (Xk, iWk) (5.34)
where fk is the system transition function and w,_, is a random, zero mean, whitenoise
sequence representing a known calculable probability distribution function. At discrete
times, an observation becomes available
z, = h(x,,vk ) (5.35)
where vk is a zero mean, whitenoise sequence representing a known distribution
function. It is further assumed that an initial pdf p(xolYo) of the state vector is available.
The objective of this process is to construct the pdf of the current state xk given all the
available information; i.e., p(xkIY).
20
The prior distribution before the update is defined as
p(xk I Y) = p(xk I x) p(xk 1 k)dxk (5.36)
where the probabilistic model of the state evolution can be obtained through the use of
the system state equations and the known statistics of wk as
P(Xk xk1) )=p(xk Ixl, wk) p(wk_xkI) dw,_ (5.37)
Because the process noise is assumed to be independent of any previous state values,
p(wk_ xk1) = p(wk_). If x,_1 and wk_ are known, then x, can be arrived at
deterministically.
The update stage occurs when a measurement becomes available. The a priori
density may be updated according to Bayes' rule as
p(xIY )= P(Yk Ix) P(xklYkI) (5.38)
P(YlY ,)
where the normalizing denominator is
P(YlYk) )= p(yYk x,) p(xk IYkY) dx (5.39)
and the conditional pdf, p(ylxk), is defined by the measurement model and the known
statistics of v,.
The method by which the following theory is implemented in the bootstrap filter
will now be defined. The bootstrap filtering algorithm provides the mechanism for
defining a set of random samples {x, (i); i = 1,..., N} from the pdf and obtaining the
distribution p(xk5, ,,_) .57 The propagation stage assumes a number of N random
samples have been taken from some initial pdf p(xk _, lY). Each random sample for
each state is then passed through the system model to obtain a set of a priori samples at
time step k
xk (i)=k 1 f (i),wkl(i)),i= 1,... N (5.40)
where the asterisk represents the a priori sample.
The update stage of the filter takes each of the a priori samples and computes a set
of a posteriori samples from which the estimate can be obtained. To do this, the
likelihood of each a priori sample l(xk'y ) is calculated using the measurement
equations and a normalized weight is given to each sample by
P(yk Xk(j))
This defines a discrete distribution over the a priori density. The normalized values of q,
can be resampled using a uniform distribution. A larger value of q, allows more samples
from that bin to fall into the a posteriori pdf. The state associated with each new sample
selected is assigned to the new a posteriori state {xk (i); i = ...,N. This a posteriori set
of samples is approximately distributed as p(xk lY). From this the mean and variance for
each state estimate can be computed. Because this filter makes no assumption of linearity
or Gaussian statistics, it was useful in providing a close approximation of the true pdf.
This semitrue pdf can then be used to investigate the effect of linearization and
parameterization on the evolution of the pdf.
22
5.3 SDREF Derivation
The SDREF is the nonlinear "dual" of a nonlinear regulator control design
technique, which uses a parameterization technique to bring the nonlinear dynamic and
measurement equations
x=f(x) (5.42)
z = h(x) (5.43)
and transform them into a linearlike statedependentcoefficient (SDC) form
x= F(x)x (5.44)
z = H(x)x (5.45)
The distinctive assumption in the SDREF technique is that the SDCs can be
treated as "frozen" over a small interval to give adequate piecewise performance. The
algebraic Ricatti equation is solved at every filter measurement update and the linear
Kalman filtering algorithms applied.
The technique was first applied to the following nonlinear regulator problem.13
The task is to minimize
I= 1 xQ(x)x + uTR(x)u dt (5.46)
with respect to the state x and control u, subject to the nonlinear differential constraint
S= f(x) + g(x)u (5.47)
The SDREF approach for obtaining a sub optimal solution is to use direct
parameterization to bring the nonlinear dynamic equation to the SDC form
i= A(x)x+B(x)u (5.48)
23
where f(x) = A(x)x and B(x) = g(x)
From this point the SDREF can be solved for P as
Ar(x)P + PA(x) PB(x)R~(x)BT(x)P+Q(x)=0 (5.49)
and the nonlinear feedback controller can be constructed as
u= R(x)BT(x)P(x)x (5.50)
The SDREF technique is shown'3 to have desirable properties of local asymptotic stability
and suboptimality given mild conditions of stabilizability and detectability.
It also shows that if A,(x) and A,(x) are distinct SDC parameterizations, they can be
combined to yield a third parameterization by
A(x, a) = a4, (x)+(l a)A2(x) (5.51)
This technique of combining parameterizations allows an infinite number of
parameterization possibilities to enhance controller design.
An example'3 of the SDC form and SDREF is repeated to introduce the method.
Consider the stochastic nonlinear system
= f(x) + r(w) (5.52)
z= h(x) + v (5.53)
where w is Gaussian zeromean white process noise with E[w(t)wT(t + r)]= Q(t)8(r)
and v is Gaussian zeromean white measurement noise with E[v(t)v'(t + )] = R(t)8(r).
After bringing the system to the SDC form
x= F(x)x+Fw (5.54)
z= H(x)x +v (5.55)
24
the SDREF is given by the steady state, linear, continuous Kalman filter equations as
x = F()i+ + K (i)[z(x) H(i)5] (5.56)
where K (j) = PHT ()R' (5.57)
and P is the positive definite solution to the SDREF
F(i)P+ PFT(i) PHT(i)R'H(i)P+F'Q = 0 (5.58)
The preceding method was applied to a simple twodimensional pendulum
problem.13 Results were comparable to those from an EKF. The main thrust of the
research was to develop a SDREF for the linear dynamic and nonlinear measurement
equations associated with the strategic, terminal homing guidance problem and
investigate filter characteristics using the EKF and bootstrap filter for comparison and
analysis.
5.3.1 AngleOnly Measurement SDREF
The following derivations are performed for the bearingmeasurements only
SDREF. It shows that this form of the SDREF is not detectable and therefore violates the
necessary conditions to assure the positive definite solution of the Ricatti equation. Both
the continuous and discrete forms of the filter are given, as well as stabilizability and
detectability analyses for each type. The filter derivation begins with the definition of the
state equations
x= Fx + Gu+ Lw (Continuous form) (5.59)
Xk =lkIX_1 +rFk_l1+Ak_Iwk_1 (Discreteform) (5.60)
and the measurement equations
z=h(x)+v (Nonlinear form) (5.61)
z = Hx+v (Linear form) (5.62)
The states to be estimated are
SXR
SYR
SZR
VXR
xk = VYR (5.63)
VZR
vz,
ayr
TaZ
az,
where the states are relative position, relative velocity, and target acceleration.
The system equations are
sR = vR (5.64)
R = AT (5.65)
a, = A + WT (5.66)
Putting this in state space notation
R = 0 L v + I aco, + O WT (5.67)
Ir 0o JAlaJ /
i = [F] Ix)+ [G] u + [L]w
The measurement equation is
Zk =h(x)+v,
(5.68)
26
where, based on the smallangle approximation,
h(x) = Rane Range = J2 +y2 R2 (5.69)
Range
and vk = N(0,R,)
The SDREF approach for obtaining a sub optimal solution is to use direct
parameterization to bring the nonlinear dynamic equation to the SDC form
x= F(x)x+G(x)u (5.70)
where f(x) = F(x)x
and z= H(x)x + v (5.71)
From this point the SDREF can be solved
F(x)P+PFT(x)PHT(x)Rl(x)H(Hx)P+LQLT(x)=0 (continuous) (5.72)
QPTgrPOPHr(R+HPHT)'HPOT+Q_, =0 (discrete). (5.73)
Since the system dynamic equations are linear, the only parameterization needed is for the
measurement equation. From the definition
z = Hx + v (5.74)
the linear form for the SDREF can be cast as
H [ R 0 0 0 0 0 0 J (5.75)
01/R 0 0 0 0 0 00
5.3.1.1 Continuous SDREF
Using the abovedefined state and measurement equations with A =0.1, the
following matrices are used for the filter equations
27
0 0 0 1 00 0 0 0
0 0 0 0 1 0 0 0 0
0 0 0 0 0 1 0 0 0
0 0 0 0 00 1 0 0
F= 0 0 0 0 0 0 0 1 0 (5.76)
0 0 0 0 00 0 0 1
0 0 0 0 0 0 .1 0 0
0 0 0 0 0 0 0 .1 0
0 0 0 0 0 0 0 0 .1
The time constant for the Markov process A, was based on an expected target
acceleration of approximately 981 (m/s2).
0 0 0
0 0 0
0 0 0
1 0 0
G= 0 1 0 (5.77)
0 0 1
0 0 0
0 0 0
0 0 0
The process noise used for this problem is 1 g on the xtarget acceleration and 3 g's on
both y and z target accelerations. The disturbance input is
E[w(t)] = 0 (5.78)
[(000)= (( )] (579)
000
000
G= 0 1 0 (5.77)
0 0 1
000
000
000
The process noise used for this problem is 1 g on the xtarget acceleration and 3 g's on
both y and z target accelerations. The disturbance input is
E[w(t)]=0 (5.78)
E[w(t),w'(V)= Qc(t)3(tr)] (5.79)
28
which for small values of At can be expressed as
0 0 0 0 00 0 0 0
0 0 0 0 00 0 0 0
0 0 0 0 00 0 0 0
0 0 0 0 00 0 0 0
Qk1 =LQcL At= 0 0 0 0 0 0 0 0 0 At (5.80)
0 0 0 0 00 0 0 0
0 0 0 0 0 0 96 0 0
0 0 0 0 0 0 0 866 0
0 0 0 0 0 0 0 0 866
where 866 (m/s2) is the variance of 3 g's of acceleration noise in the y and z directions.
The parameterized measurement matrix can be computed as
0 0 0 0 0 0 0 0
(x)= Range (5.81)
0 a 0 000000
Range
The standard linear, continuous Kalman filtering equations follow as
= F(i) + K, (G)[z() H(i)i] (5.82)
where K, (i) = PHTR,' (5.83)
and P is the positivedefinite solution to the statedependent Ricatti equation
F(i)P + PFT() PHT(i)RK,_1H(i)P+Q, = 0 (5.84)
where Qki = LQc'LrAt. Ricatti equation theory provides for a legitimate positive semi
definite solution if the following conditions are met
(TD,Q,PER",HT ERm,ReR RrnQ=QT>0,R=RT>o,m
and (F,H) is pointwise detectable in the linear sense and(C, F) is pointwise
29
stabilizable in the linear sense for all x e 9, where C is a fullrank factorization (FRF)
of Qk given by
(C'C)= Q and rank(C)= rank(Qk) (5.86)
For the given linear dynamics we can construct the controllability matrix as
Con=[C FC ...F"'C] (5.87)
If Con is of full rank, controllability and therefore stabilizability is assured. If it is not
fully controllable, a minimum condition of stabilizability must be met. Putting the
system equations into the controllability canonical form, the uncontrollable poles can be
checked for stability.20
Using a similarity transformation, the controllability canonical form can be found
as
F=TFT', C =TC (5.88)
where T is unitary and the transformed system has a staircase form
= 0] C_ 0o
where (F,,G.) are controllable and (F,,) is uncontrollable. The poles of the
uncontrollable subspaces are then checked to insure stability. If they are found to be
stable, then stabilizability is verified.
Using the noise matrix in equation (5.13) and performing the FRF for C, it can be
shown that the controllability matrix has full rank and therefore stabilizability is verified.
30
Similarly, the observability Grammian matrix is constructed as
H
HF
Obs= (5.90)
HFn1
If Obs is of full rank for all x, then full pointwise observability in then linear sense is
assumed. Once again pointwise detectability needs to be assumed for the Ricatti
equation solution to be positive semi definite. If it is not of full rank the similar
observability canonical form can be used to again check the pointwise stability of the
unobservable poles. This canonical form is given by
=[F FI, H=[0 Ho] (5.91)
Using the above referenced matrices for Fand H, it can be shown that the observability
matrix is not of full rank. Therefore, pointwise stability of the unobservable poles must
be checked. It was found that the poles were not pointwise stable and therefore point
wise detectability was not verified. A similar analysis was then performed on the discrete
form of the filter. From this point on, only the discrete form of the SDREF is analyzed.
This form was chosen due to the use of discrete measurements that are only available at a
subset of the integration step size of the simulation. The discrete filter is better suited for
this type of simulation and would most likely be used in any realistic guidance software.
5.3.1.2 Discrete SDREF
The discrete form of the SDREF is derived using the following formulas. Since
F is a constant matrix, the state transition matrix can be found using Laplace transforms
as
I IAt J(e +At1)/1
D(At) = 0 I I(1e ')/A (5.92)
0 0 le
When using 2 = 0.1 and At = 0.025, the following transition matrix is computed
1 0 0 .025 0 0 .000312 0 0
0 1 0 0 .025 0 0 .000312 0
0 0 1 0 0 .025 0 0 .000312
0 0 0 1 0 0 .0249 0 0
S= 0 0 0 0 1 0 0 .0249 0 (5.93)
0 0 0 0 0 1 0 0 .0249
0 0 0 0 0 0 .9975 0 0
0 0 0 0 0 0 0 .9975 0
0 0 0 0 0 0 0 0 .9975
Assuming the commanded acceleration remains constant over the filtering interval
At2A co
2
Fuk ,,= (tk,T )G()ud= AtAcoml (5.94)
0
The state propagation equation becomes
[ At21
2
x [E]x, + AtI Ao (5.95)
0
Also,
E[wwT] = QS(t ) (5.96)
and
Qk, = I(tk ,T)L()Q'c (r)L (r)c'(tk ,r)dv (5.97)
1I
Using 1 g of acceleration noise on the x target acceleration and 3 g's in y and z and
performing the above integration, the process noise matrix is obtained.
9.46e9 0 0 9.36e7 0 0 4.99e5 0 0
0 8.54e8 0 0 8.44e6 0 0 4.45e4 0
0 0 854e8 0 0 8.44e6 0 0 4.45e4
9.36e7 0 0 9.98e5 0 0 .006 0 0
Qk = 0 8.44e6 0 0 9.0e4 0 0 .054 0
0 0 8.44e6 0 0 9.0e4 0 0 .054
4.98e5 0 0 .006 0 0 .478 0 0
0 4.5e4 0 0 .054 0 0 4.32 0
0 0 4.5e4 0 0 .054 0 0 4.32
The discrete Ricatti equation solver requires the following conditions to be met to insure
that a legitimate positive semi definite solution to the Ricatti equation is obtained:
(O, Q,Pe Rn',HT R Re Rm", Q=Q > R = R >0,m
and (D, H) is pointwise detectable in the linear sense, and (C, () is pointwise
stabilizable in the linear sense for all x e Q, where C is a FRF ofQ,_ and is defined as
(CTC)= Qk and rank(C)=rank(Qk_) (5.99)
It is also must assumed that Q is invertible without need for further modifications to the
solver. If all assumptions are met, then the Ricatti equation will have a unique non
negative definite solution throughout 2.
33
Using the above referenced matrices, the same analysis used for pointwise
stabilizability of the continuous system was performed. The controllability matrix of the
pair (C,D) was found to be of full rank and therefore the system is stabilizable.
Again an analysis similar to the continuous version was performed using the pair
((D, H) and the pointwise observability matrix was again found to be not of full rank.
The canonical form was calculated and again the poles of the unobservable subspace were
determined to be pointwise unstable and therefore the system could not be declared
pointwise detectable in the linear sense.
Because the system has been shown to be pointwise undetectable, the Ricatti
equation solver cannot provide the positive semi definite solution to the Ricatti equation.
Consequently, another approach to the problem must be examined. One way to make the
system detectable is to allow a measurement in the range or xdirection. This is shown in
the next section.
5.3.2 SDREF With Range Measurement
This section analyzes the same filter described above except that a range
measurement has been added to the measurement equation. By adding the range
measurement the system becomes both stabilizable and pointwise detectable under the
conditions set forth above. Two types of range measurements were evaluated to provide
adequate observability for the system.
5.3.2.1 Intensity range measurement
The range measurement used is the source radiant intensity of the target arriving
at the focal plane and is
I= Range k (5.100)
Range
where
Ik =intensity (watts)
( watts i
J, = source radiant intensity watts
(sr/pm)
A = bandwidth (,um)
r = radius of seeker aperature (m)
R = range from source to seeker (m)
v, = white Gaussian noise
This gives the measurement equation:
z/Range
h(x)= ylRange Range=x x/+yR2 +zR (5.101)
JA2tmr2 /Range2
Putting equation (5.101) into the SDC form
0 0 000000
Range
H= 0 1 0 0 0 0 0 0 0 (5.102)
Range
Rang 2r 0 0 0 0 0 0 00
Range2 *XR
The measurement noise for the range measurement is somewhat more complicated as it is
range dependent. The actual intensity measurement has Gaussian noise from the focal
plane electronics, which can cause error in the measurement. However, it is assumed this
error is small in comparison to the fluctuations of the source radiant intensity and is
35
assumed to be zero. Jz is, however, comprised of both the deterministic source radiant
intensity and a random part given by
JA = J, + I_ (5.103)
The resulting intensity contribution from the random fluctuations in the plume
core can be calculated from equation (5.101) as
I = J (5.104)
Range2
The variance of this random component of intensity is obtained by taking the expected
value of the square of I, (i.e., E{I, }). The equation for the expectation operation
gives
E i(J,,AAIr2' 2 a E{J' 1 (5.105)
Range1 J Range2
If p= E(J2,, ), then the range dependent measurement variance term for the filter
becomes
(pAArr2 2
SRange 2 (5.106)
Using appropriate values for these constants, a measurement matrix can be given by
2.5e9 0 0
Rk= 0 2.5e9 0 (5.107)
0 0 8.9e 7
Using the abovedetermined matrices, the SDREF was again calculated. The filter
responded well during the first three to four seconds of flight and then diverged. The
36
problem lies in the instability of the Ricatti equation solution. In examining the point
wise closedloop eigenvalues of the system (i.e., the spectrum of
DT HT(R + HPHT) HPQ (5.108)
with x, in this case range, frozen), the system would be asymptotically pointwise stable
in the linear sense if the closedloop eigenvalues all lie inside the unit circle. The
eigenvalues of the SDREF were computed for each filter update stage. They began near
the unit circle (i.e., 0.9992) and grew to become unstable as the simulation ran. The final
step in the Ricatti equation solution is the solution of an nt order linear matrix equation.
A number was computed to determine the condition number for the inverse matrix. It
could be seen that the condition number increased (i.e., illconditioning) as the closed
loop eigenvalues approached the value of one. The relatively small size of the range
variance causes R' to become very large, and therefore ill conditioned. Once the closed
loop eigenvalues became larger than one, the problem diverged quickly. Because of this
numerical instability, a different range measurement was tried.
5.3.2.2 Distance range measurement
The new range measurement is the distance measurement in meters. This value is
assumed to be available from the signal processing algorithms. The new measurement is
Range =x, + yR + (5.109)
37
Therefore, the new measurement equation becomes
z
Range
h(x) (5.110)
Range
Range
and the measurement matrix becomes
0 0 1 0 0 0 0 0
Range
H= 0 0 0 0 0 0 0 (5.111)
Range
Range 0 0 0 0 0 0 0 0
X,
Using these values in the SDREF, satisfactory results were obtained. (These results are
given in Chapter 7.) It was decided that this would be the formulation that would be used
to further test this new filtering technique. Before these results are presented, however,
an investigation into the theoretical stability of the SDREF is presented.
5.3.3 Theoretical Investigations
One of the main theoretical criteria for filtering techniques is the stability of the
algorithm. The following section reviews several theoretical stability investigations.
While these investigations were unsuccessful in proving stability, it is hoped that they
provide valuable insight into the SDREF. It is these attempts that led to the development
of the MSDREF shown Chapter 9. While the methods outlined below approach the
problem from several different theoretical aspects, all were unsuccessful in showing local
stability for the SDREF. Because the true and estimated states x and i both appear in the
error equation
38
e = (F(x)x F(.)! K(x)[H(x)x H(i)]) (5.112)
each of the approaches fail when presented with the cross terms that are created. It is
hoped that this section will prove useful to future stability investigations. The approaches
are shown starting with simple linearization techniques and progress through more
complicated and sometimes restrictive theorems. The first of these is a simple Taylor
series expansion of the parameterized equations.
5.3.3.1 Taylor series expansion
The following stability analysis will exploit a Taylor series expansion about the
local zero point. Consider the stochastic, nonlinear system
= f(x)+g(w) (5.113)
z= h(x)+v (5.114)
where w is Gaussian zeromean white process noise with E[w(t)wr(t + r)] = Q(t)8(T)
and v is Gaussian zeromean white measurement noise with E[v(t)v'(t + )] = R(t)8(T).
After bringing the system to the SDC form
= F(x)x+Fw (5.115)
z = H(x)x+v (5.116)
the SDREF is given by the linear, continuous Kalman filter equations having state
dependent coefficients
S= F(i) + K(i)[z(x) H(i)i] (5.117)
where K(i)= P(i)HT(i)R' (5.118)
and P is the positive semi definite solution to the algebraic SDRE
F(i)P(i)+P(i)FT(i)P(i)HT(i)R'H(i) P(i)+ FQF =0 (5.119)
39
Conjecture: Assume that R > 0, Q >0, f(x), and h(x) RK, and that the
parameterizations F(x) and H(x)can be shown to satisfy that (F(x), H(x))is a pointwise
detectable pair in the linear sense and (C(x), F(x)) is a pointwise stabilizable pair in the
linear sense, where CC = Qk,. Also assume F(x) e C" and H(x) e C~ are analytic,
where C~ is the class of functions whose derivatives of all orders are continuous, and that
f(0) = h(0) = 0. Then the SDREF is locally asymptotically stable.
Outline of Attempted Proof:
For ease of notation the proof will be examined in scalar form. Given that the
scalar system dynamic equations are
x = F(x)x+Gu+ Lw
z=H(x)x+v and (5.120)
x = F(i)i + Gu +K[z H(x) ]
and defining the estimation error rate as e = x the above equations can be written as
= (F(x)x F(i) K(x)[H(x)x H(i)i]) (5.121)
Now expanding h(x) and H(x)x in a Taylor series expansion about zero
h(x) = H(x)x (5.122)
A(0) d h ) h (0) i (x0)
(0) (0) + H.O.T H(0)+ )+ +.O.T x
x= 10 A X=o 2! L Ix=0 ax=0 21
Also for f(x)= F(x)x
f(x) F(x)x (5.123)
3() d2f ( 2O)2 F dF(O) d2F (x0)2
Po)+ (xo)+ if + 0) o L H O.T
(0o )+ +H.o.T.= F(O)+ ((x)+" 2 o.r. x
SI o= 2 L A X=O a 0 2!
Similar equations can be written for h(I) = H().i and f (x) = F(i)i.
Equating like terms
A(0) A(0)
h() = H(0) 
and O(o) F(O)= (O)
Using this and recalling that K = PHTR',
E(t) = (F(0) PHT(O)R'H(O)))E(t)+ H.O.T. (5.124)
Herein lies the problem, the higher order terms contain both x and i. These terms cannot
be guaranteed to approach zero, much less go to zero faster than the first term. If the
error dot equation could be expressed in this linear form a local stability proof could be
achieved in the following manner. For ease of notation, let F(0) = F, H(0) = H, and
e(t) = e. Following a similar derivation,'9 a Lyapunov function is chosen
V(e) = eIE (5.125)
where I is the positive definite information matrix, the inverse of P. From the definition
of the derivative of the inverse,
I = P = IPI (5.126)
By substituting this into the matrix Riccati equation, a differential equation
j = IF + FI +ILWLTI HTN'H (5.127)
is yielded that has the steadystate solution
0= IF +FrI +ILWLI HTN'H (5.128)
41
The time derivative of the Lyapunov function is
V= 2erI = 2eT[F KH]e
= 2eT[F PHTN'H]e (5.129)
= er[(IFHTN'H)+(IF HTNH)]e
Using (5.128) this can be written as
V= 4ILWLrI+H TN H]E (5.130)
As the term in square brackets is positive definite at all times, V < 0 at all times. This
would verify the local asymptotic stability of the SDREF for some small neighborhood
about the origin. The estimation error would approach zero as time approaches infinity.
5.3.3.2 Using a Taylor series expansion about the true states
This approach varies slightly from the above, in that a Taylor series expansion of
the estimated dynamics and measurement equations are linearized about the true states.
Using the same state and measurement equations (5.1135.120) and the error equation
(5.121) in the preceding proof, F(i) and H(i) are expanded as
F"(x)
F() = F(x)+F'(x)(ix) + (.x)2+...H.O.T. (5.131)
2
and H()= H(x)+H'(x)(x)+H(x x)2+...H.O.T. (5.132)
2
Using equations (5.131) and (5.132) equation (5.121) is written as
S= {F(x) KH(x)}e+ {F'(x) KH'(x)} d t(e2) (5.133)
Once again the problem with this approach can be seen. The only way it can be shown
that the second term disappears or is negative is for ex to be zero. This can only be
shown to be true for the linear case.
42
5.3.3.3 Taylor series expansion about the error
Again, using the same system equations (5.1135.121) as in the first Taylor series
approach, this second approach linearizes about the error. Linearizing the dynamic and
measurement equations
F"(E)
F(x)= F(e)+F'(e)(x )+ (xE)2+...H.O.T. (5.134)
2
F"'(e)
F'(x)= F'(e)+F"(E)(x e) + (x )2+...H.O.T. (5.135)
2
H"(e)
H(x)=H(e) +H'(E)(x ) + (x e) +...H.O.T. (5.136)
2
H"(e)
H'(x)= H'(e)+ H"(e)(xe)+ (xe)2+... HO.T. (5.137)
2
After substituting these into error equation (5.121),
= (F KH) +(F' KH')(x e)e
+(F' KH')ei +(F" KH")(x E),
+(F"KH")(x E)2 (5.138)
2
+(F"'KH"')(x e)2 + (e2)
2
Again, the error equation is only stable if it can be shown that ex > 0. This can only be
shown to be true for the linear case.
5.3.3.4 Residual error stability approach
The following approach illustrates, under conditions of pointwise detectability in
the linear sense that ensure a legitimate solution of the SDRE, that the SDREF can be
shown to be asymptotically stable.21 The system dynamics and filter equations are given
in discrete form as
xk+i = 4kXk +Wk (5.139)
Zkt+ = h(xkl)+ Vk+ (5.140)
Xk+i/k b (5.141)
"k+i =k+k + Kk+I[Zk+ Hk+ (k+Ik)k+l/ ] (5.142)
Kk+I Pk+1 (k+I)Hk+T (k+,)Rk+,' (5.143)
with E[vk+l vk+I ]= Rk+l (5.144)
E[wwk+WT] = k+ (5.145)
where P, (ik+) is the solution of the SDRE. These are the equations necessary for the
discrete SDREF. It is also assumed that vk+5 and w.k+ are random, zero mean inputs and
therefore, the stability of the SDREF is independent of these error sources.
The state estimation error both after and before the update is defined respectively
as
k k+ k+ (5.146)
xk+/k = Xk+l Xk+/k (5.147)
Also, for ease of notation
Hk+ (k+= Hk (5.148)
Pk+ (k= P+. (5.149)
e,+ = Zk+l Hk+ (k+ )k+l (5.150)
And finally a candidate Lyapunov function Vk, is defined as
vi = pk+ pk Xk+ (5.151)
The conditions for which Vk+~ is a decreasing sequence must be determined, therefore
showing the exponential stability of the SDREF.
44
In the preceding stability approaches, linearization of the filter equations has not
shown stability of the filter. This attempt constrains the state transition matrix and
measurement noise matrix so as to characterize the error due to linearization. The
measurement residual error can be characterized as
ek+1 =zk+ h(ik+1/k) (5.152)
e+ = Hk+(l k+ )Xk+l Hk+ (k+ l)k+l (5.153)
Rewriting H,, (i+,) using a Taylor series expansion, for each output error prediction
component of e+,, namely e, ,k+, there always exists a residual r,,k+jdue to the first order
linearization in order to obtain an exact equality. Therefore,
e,.k+1 H,,k (x+ k+I)x+ Hi.k+ (xk+) + H'i,+ (xk+1 )(X+ X )+. (5.154)
and the residual error is
rk+= {Hi, (Xk+1)(ikl Xk)+. } +l (5.155)
Then the exact measurement error can be defined as
e,.k+ = Hik+1 +1k + r.,k+ (5.156)
or written in another form
Hik+3lk+1k = ai.k+lei,,k+ (5.157)
where a,,,k is an unknown timevarying scalar related to ri,k+ by
,k++ (1 a,k+l )e,k+l (5.158)
Introducing the signal vector
ak+lek+ = Hk+15k+/k (5.159)
where ak+, e RP' is an unknown time varying diagonal matrix defined as
ak+1 = diag(a,k+1, .. a.,~k+1) (5.160)
If both sides of equation (5.142) are subtracted from xkl,,
+ = k+I/k Pk+IHk+ITRk+I'ek+1, (5.161)
Substituting this into equation (5.151), and using the fact that Pk+1 =Pk
Vk+1= .k+k iT Pk+11 k+l/k Ktk T Hk+RIT+ k+1 (51
ek +lt R+1 Hk+1 k+l/tk ++l tR+ Hk+,Pk+H l R,+ ek+l
Because the covariance matrix is not propagated during the time between updates, it is
assumed that the a priori covariance before update is equal to the previous updated
covariance,
Pk+,1k = P (5.163)
Therefore,
Pk+' =Pk' (5.164)
Before proceeding further, the inverse relationship is defined as
Pk+ = Pk+,,11 + Hk+,TRk+,'Hk+5 (5.165)
Once again a similar problem arises. Equation (5.165) is not derivable from the basic
premise of the SDREF update equations.21 If this equation was available, it could be
shown that Vk,, is a decreasing sequence. Then
Vk+I = Vk+I1, Xk+U III T H Rk+1 ek+1
(5.166)
ek+, Rk+1, Hk+,,k+11k +e,, +RLHk+ Pk+ Hk+IT Rk+1 ek+
and reviewing that
Hk+lk+1= ak+, (5.167)
and
XTk+IHHTk+i = eTk+laT k+
(5.168)
Equation (5.166) is written as
Vt+l = Vk+l/k ek+,1 t+a k+l atk+lek
ek+ ak+lRk+lHk+lek+1/k k e+TRk1 'ak+ek+ (5.169)
+ek+l Rk+l lHk+I Pk+lHk+I Rk+,lek+l
And using the fact that
Vk+,1/ = TFkTPk+,,kI'Fk (5.170)
with Equations (5.167168), Equation (5.169) is rewritten
Vk+1 V = eV+lT k a+l Rk+,1ak+l akt+iRk+., Rk+,,ak+,
+Rk+1lHH+lP+1Htk+TRk+'I }ek+, (5.170)
+;5FT(FkTPIF Pk1
Therefore under the condition of pointwise detectability in the linear sense, and certain
restrictions on R and Q (i.e., Rk > 0 and Qk >0) then,
VI, <0
and the system is asymptotically stable as the values in parentheses are < 0.
5.3.3.5 Lyapunov approach using the expectation operator
This stability approach uses a similar method.9 Given the system equations
x+, = A, ,+w, (5.171)
= + K,[H(x,)x, H(x,)x,] + Kv, (5.172)
the error can be updated and propagated by the equations
e, x, i (Updated error) (5.173)
x x, (Propagated error) (5.174)
Using equations (5.1715.172)
e, = x, x, K,[H(x,)x, H(,),] K,v (5.175)
= A,_e, + w,_, (5.176)
After some manipulation of these equations, the error update equation can be developed
eF1= A, , Kix,+A,K,HA, H K,v,+w, (5.177)
where H, = H(x) and H, = H(()
If the following Lyapunov function is defined
V,+ = e,j.e+, (5.178)
it can be shown that
E{V, (eF.)V,(e,)} 0 (5.179)
Substituting equation (5.177) into equation (5.179), the following equation can be
computed
E{V,,,((,) V()}= eE{ATA I}e,
+E{2"ATAKH,x2xTATAKHx}
+E{.2xTArAKH2x + 2xTATAKH2Y}
T T KT AT Y1 (5.180)
+E(2xTHTKTATAKH (5.180)
+ ExTHTKTATAKHlx}
+ E{xTH2,TKTAAKH2x}
In this equation several terms exist that cannot be shown to be negative definite. These
terms are brought about by the same problem found in the previous approaches evaluated
to show stability. They are cross products created by the presence of both xand i in the
error equation.
This same approach was also tried with the Lyupanov function
V=eTe (5.181)
and Equation (5.182) was developed.
E{V,, (e,)V(e,)} = e,TE{ATA e,
+ E{2A'rAKH,x 2x ATAKHx}
+ E(2 TArAKH2i + 2xTA'AKH2i}
+E{2xTHTKTATAKH2} (5.182)
+E{XTH,TKTKKH,x}
+ E(IrH2T'KKH2z)
Once again cross terms are generated and cannot be shown to be negative definite. All of
the previous attempts at stability proofs fail for the same reason: The cross term xi,
produced in the error equation, raises significant difficulties in finding a closed form
differential equation for e. Attempts to limit the problem to a scalar were not able to
resolve the cross coupling problem. While several standard and nonstandard approaches
were tried, a satisfactory stability proof could not be verified. It is hoped that these
investigations can prove useful in future attempts to determine filter stability.
CHAPTER 6
EXOATMOSPHERICGUIDANCEPROBLEM SIMULATION RESULTS
Chapter 6 investigates the comparative performance of each filter against a variety
of engagement conditions and noise parameters. Before any study on filter performance
was done, a baseline performance analysis was completed, which allowed a detailed
investigation of the SDREF against one engagement scenario. Once the baseline
performance was established, further investigation into filter sensitivity, robustness, and
closedloop performance was performed. A study of pdf evolution for each filter is also
shown.
6.1 Baseline Performance
To establish a baseline performance measure, an engagement representing the
final 10 seconds of a closing interceptor/boosting ICBM engagement was chosen. Table
1 details the major initial filter parameters. The aspect angle is defined as the initial Euler
rotations between the axes of the interceptor and the target. For example, an aspect angle
of (30,30,30) would signify that the interceptor's initial velocity vector has been rotated
30 degrees about each axis. An aspect angle of (0,0,0) would be a headon engagement.
The (30,30,30) aspect angle was chosen so that all axes of the filter would be exercised in
the engagement. SIGP and SIGVP are the standard deviations for the random
initialization errors on the "true" position and velocity states respectively. When SIGP
50
and SIGVP are zero, the filter has been initialized with the "true" initialintercept
conditions. SIGXHAT are the values used to initialize the covariance matrix for the EKF
and the bootstrap filter. XHAT is the initial state given to the filters. The measurement
noise for the xdirection is a rangedependent noise that is given as a percentage of the
range (i.e., as the interceptor gets closer to the target the noise on the range measurement
decreases). All filters were made to operate with the same initial conditions and statistics
so that comparisons could be made.
A perfect filter was used to guide the vehicle to impact, while the EKF, SDREF,
and bootstrap filter were run open loop. The bootstrap filter was run with 10,000
samples.
Table 1. Filter simulation parameters baseline conditions
Parameters (units) X Y Z
Aspect Angle (degrees) 30 30 30
SIGP (m) 0 0 0
SIGVP (m/s) 0 0 0
a,2 (m2) 1.0E5 1.0E5 1.0E5
SIGXHAT O 2 (m2/s2) 1.0E3 1.0E3 1.0E3
aA_2 (m2/s4) 1.0E2 1.0E2 1.0E2
R (m) 1.0E5 10.0 10.0
XHAT V (m/s) 1.0E4 335.9 223.9
A (m/s2) 55.2 57.40 38.2
SIGQ (EKF) (m2 /s4) 196.0 196.0 196.0
SIGQ (SDREF) (m2 /s4) 1.0E6 1.0E3 1.0E3
Meas Noise (%,pm,/pm) 0.1 50.0 50.0
51
The following graphs reflect filter performance. Figures 3 through 5 show the
positionstate estimate errors with 3v+. standard deviations for each filter. Velocity and
acceleration performance was similar between the EKF and SDREF. The bootstrap
acceleration estimates are poor due to an inability to forward measurement information to
the second derivative states. The problem can be remedied by the use of millions of
samples or other adhoc improvement schemes; however, computer limitations and the
scope of this dissertation do not allow investigation at this time. Figures 3 through 5 give
initial insight into the relative performance of each of the filters for the given scenario.
All of the filters adequately estimate the state to the precision necessary for an intercept to
occur. With state estimate errors remaining within 3a standard deviations, the filters
perform well. The bootstrap filter does show a tendency to keep a wider variance than
the other two filters. However, as discussed in paragraph 6.2, when the EKF and SDREF
are run in Monte Carlo fashion, they too have wider variances. The results shown here
are for one run, while noticing that the bootstrap filter is inherently Monte Carlo, even for
one run. The performance curves in Figures 3 through 5 verify proper filter operation and
establish a baseline performance from which more indepth analysis proceeds.
It is interesting to note in Figures 6 and 7 that the standard deviations for the EKF
and SDREF are approximately identical after the initial transients in the EKF die out.
This indicates that once the EKF has settled and the estimates are being computed
accurately, the EKF covariance propagation and update are equivalent to the algebraic
Ricatti solution being computed in the SDREF. This relationship is only shown for one
state; however, similar results were found for all states.
Figure 3. EKF position errors with +3" standard deviations.
S. 1n
,~~0~~~ I I I
=I
Figure 4. SDREF position errors with 3ustandard deviations.
1 & I I I I
'*
E
n ~m
Pm
i
0.
Figure 5. Bootstrap position errors with 3astandard deviations.
'
"
I <
!
'
r
k
10
 
Figure 6. Comparison of EKF and SDREF covariance
3cr standard deviations for Y position.
5 EKF
EKF
1 SDREF
SSDREF
15 _____ ___
15 0. 0 1
5 0.0 0.5 1.0
Time (sec)
Figure 7. Expanded view of Figure 6.
1.5 20
54
6.2 PDF Distribution Comparison
Using the baseline scenario, a comparison of the posterior (after the measurement
update) probability density functions was investigated, to determine whether
parameterizing versus linearizing the nonlinearities has an effect on the shape of the pdf.
The scenario given in Table 1 was run once and the state estimate and standard deviation
were computed three times during the engagement. Samples were taken at 1, 5, and 9
seconds. The pdfs for the EKF and SDREF were computed using a Gaussian
distribution function based on the mean and standard deviation at each time increment.
The resulting curves, Figures 8 through 9, do not represent the "true" pdf as nonlinear
dynamics are involved. However, they do show the basis of the theoretical filter
formulation. ("True" pdf functions, calculated using Monte Carlo simulations, are
discussed in following paragraphs.)
The pdf for the bootstrap filter was generated by computing a histogram of 1,000
samples, curve fitting a function using a cubic spline, and normalizing the area under the
curve. Figures 8 through 9 show the onerun posterior pdf functions for the filters. The
first time increment is taken 1 second into the flight. The filters are still trying to
converge and therefore the distributions are large. This is demonstrated in Figure 8 by
comparing the narrow spike of the highly converged bootstrap filter to the wide, non
converged pdfs for the other filters. If consideration is given to the scaling on the xaxis,
no concern should be given to the fact that the distributions are not near the true state. It
can be seen that the estimates are close to the "truth". The xposition state shows a more
accurate representation of the pdf than the y and zposition states. This is due, in part, to
measurement noise decreasing as the engagement progresses, and a more accurate
True Sta p
0.025  oostrap
Time= 1 Second EKF
SDRE
0.020
0.015
S0.010
0005
0.000
88500 89000 89500 90000 90500 91000
X Position (m)
0.05
Time =5 Seconds True sate
0,04
0.03
S0.02
0.01
0.00 ... . ...
50000 50200 50400 50600 50800 5
X Position (m)
0.16
0'14 Time 9 Seconds True St 
 EKF
0.12
0,10
t 0.08
008
0 0.06
0.04
1000
12180 12200 12220 12240 12260 12280
X Position (m)
Figure 8. Posterior probability density functions at
time = 1, 5, and 9 seconds for X position.
I
Y Position (m)
True State
5 T t Boolstrap
S EKF
Time = 5 Seconds SDREF
1
421 420 419 418 417 416 415
Y Position (m)
 Bootstrap
EKF
Time = 9 Seconds SDREF
True State
174.0
1735 173.0
Y Position (m)
Figure 9. Posterior probability density functions at
time = 1, 5, and 9 seconds for Y position.
0.
0.
0.0
5
4
3
2
knowledge of the true state is obtained. Another interesting factor shown in Figure 9, is
the compression of the EKF and SDREF distributions. This is typical of the "falling
asleep" problem associated with many filters. As the state estimation error becomes
small, the covariance matrix also becomes small, and the filter reacts slowly to new data,
which takes the estimate away from the converged answer. The bootstrap filter, however,
keeps the distribution open wider. Even using the onerun pdfs, the EKF and SDREF
pdf shapes appear to be similar.
To determine the "true" pdf functions for the EKF and SDREF, the filters were
run in a 100run Monte Carlo experiment. Filter state estimates were saved at the same
time increments. The samples from each of the runs were used to generate the pdf curves
in Figures 10 through 11. The same histogram method described in the previous
paragraph for the bootstrap filter was used to plot the curves. (The bootstrap filter did not
need to be examined in this fashion, as its formulation is inherently Monte Carlo.) The
largest differences can be seen in Figure 11. In the later time samples (5 and 9 seconds),
the previously tight distributions have flattened out, while the wide distributions have
tightened up. This is to be expected, as the onerun case is only one realization of the
Monte Carlo runs.
Summarizing the above analysis shows that the EKF and SDREF pdf's do not
resemble the "true" pdf. However, they are similar to each other. This may be explained
by the fact that the EKF linearization method provides an accurate filtering methodology
as long as the filter state estimate is close to that of the "true" state. If a target maneuver
is introduced that causes the state estimate to be far from the "truth", then it is assumed
0 2798 89740 897 0
89780 89800
X Position (m)
89820 89840 89860
X Position (m)
12220 12240 12260
X Position (m)
Figure 10. Monte Carlo posterior probability density functions at
time = 1, 5, and 9 seconds for X position.
True Sale BoosIrap
EKF
Time = 1 Second SDRE
,4 '
0.05
0.04
0.03
0.02
0.01
0.00
0.35
0.30 Trne State oottp
Time= 1 Second SORE
0.25
020
0.15
0.10
0.05 
0.00
270 260
Y Position (m)
0.5
True Slate
0.4 Time = 5 Seconds Sa
03
S0.2
0.1
0 0 .. .
425 420 415 410 40
Y Position (m)
Bootsrap
EKF
Time = 9 Seconds SDREF
True State
0
175.0 174.0 173.0 1720
Y Position (m)
171,0 170.0
Figure 11. Monte Carlo posterior probability density functions at
time = 1, 5, and 9 seconds for Y position.
60
that the distributions between the EKF and SDREF differ due to the different theories of
linearization versus parameterization. The following paragraphs explore this scenario.
The filters were once again run in a MonteCarlo simulation. However, this time
the target performs a thrust cutoff maneuver, 2 seconds after launch, which causes the
target acceleration, in all axes, to drop to zero for the rest of the flight. This is
representative of timing the engagement too close to booster burnout. The filters are still
operating under the assumption of a GaussMarkov target acceleration model, even
though no target acceleration is present after 2 seconds. The same initial conditions listed
in Table 1 were repeated.
Figures 12 through 14 show the divergence from the "true" accelerationstate for
the filters. For previously discussed reasons, Figure 14 shows the bootstrap filter
performs poorly when tracking the acceleration states. The deviation from the "true"
state is significant after the maneuver occurs. This is a good test to determine the
differences between linearization and parameterization. Once again the Monte Carlo pdf
functions were plotted in Figures 15 through 16 Samples were taken at 2, 2.5, and 3
seconds into the engagement. This allows the pdf to be examined before and during the
most stressful time of the maneuver. Suprisingly, the EKF and SDREF pdf shapes do not
differ significantly. It can be seen in Figure15 that although the "true" probability has
become bimodal, the shapes are still similar. While the singlerun EKF and SDREF
distributions can not be bimodal, Figures 15 through 16 represent Monte Carlo results.
The results from the pdf comparison experiments indicate that there is no
appreciable difference in pdf shapes between the EKF and SDREF. The conclusion may
61
 I I
Figure 12. EKF estimated and true acceleration.
Figure 13. SDREF estimated and true acceleration.
Figure 14. Bootstrap estimated and true acceleration.
0 0997 79920 799 0
74900 75000 75100
X Position (m)
75200 75300
70100 70200
X Position (m)
Figure 15. Monte Carlo posterior probability density functions at
time = 2, 2.5, and 3 seconds for thrust cutoff X position.
True Sta
Time = 2 Seconds / :'
.i

0.030
0.025
0.020
N 0.015
0.010
0.005
0U000
I8ZU 79 /840 79860 79880
X Position (m)
Time= 2.5 Seconds  ooltrap
True Bate. E
\
74800

0.35
0.30
0.25
0.20
, 0.15
oto
010 
005
0.00
390
Tue
Time = 2 Seconds
 ootstrap
EKF
REF
Y Position (m)
405
400
Y Position (m)
0.40
Bootstrap
0.35 Time =3 Seconds EKF
SDREF
0.30
0.25
True State
0.20
S0.15
0.10 /
0.05
0.00 
404 402 400 398 396 394
Y Position (m)
Figure 16. Monte Carlo posterior probability density functions at
time = 2, 2.5, and 3 seconds for thrust cutoff Y position.
be drawn that the two filters operate with identical statistical models although the filters'
method of handling and operating on the stochastic processes are significantly different
6.3 Sensitivity Analysis
The following section examines the SDREF sensitivity to measurement noise,
initialization error, and the robustness properties in the presence of stressing target
maneuvers. Only the EKF and SDREF are used for this area of the study. Further
bootstrap filter development is needed to make it useful in less than ideal conditions. The
EKF and SDREF were tuned separately to give them the best chance of performing well
under these varying conditions. Therefore each filter may have been started with different
initial conditions and statistical information. The following studies are accomplished
using perfect state information to guide the vehicle, while the EKF and SDREF run open
loop. Allowing the guidance algorithms access to perfect state information removes the
closedloop guidance algorithm effects from the analysis of filter performance.
6.3.1 Initialization Error
To test the response to initialization errors, the filters were given erroneous initial
state estimate errors of 1 and 5 percent. The EKF initial covariance was set to
approximate this error. The SDREF, however, has no direct means of initializing the
covariance matrix as it is solving for the covariance matrix at each filter update time.
This will be shown to be a detriment for this particular implementation of the SDREF.
Table 2 gives the initial filter conditions for the EKF and 1 percent initial error.
Measurement noise was held constant for both cases, as was the measurement noise
matrix given to both filters. Range measurement noise is again range dependent and is
65
given as a percentage of "true" range. The initial setup for the SDREF used the same
parameters as those of the EKF except that higher values of process noise were tried to
increase the covariance response. The only tuning parameters in the SDREF are Qk, R,
and the nonlinear parameterization methodology. The measurement noise matrix,R, is
not directly related to the initialization error; therefore, only the other two tuning
parameters can be used. Because this implementation of the SDREF has linear dynamic
equations the only method available for tuning is to adjust Q,,. This is an unacceptable
way of tuning the filter for initialization errors as it causes the covariance estimate to
remain artificially high for the entire flight. The SDREF as previously derived, does not
permit any initialization of the covariance matrix, as it is simply calculated from the
solution of the SDRE. This deficiency is examined later in this section.
Table 2. EKF and SDREF parameters onepercent initialization error
Parameters (units) X Y Z
Aspect Angle (degrees) 30 30 30
SIGP (m) 1000 35 30
SIGVP (m/s) 100 3.4 2.2
o72 (m2) 1.0E6 1.0E3 1.0E3
SIGXHAT 0,2 (m2 /s2) 1.0E4 50.0 50.0
oA2 (m2/s4) 1.0E1 1.0E1 1.0E1
R (m) 1.0E5 10.0 10.0
XHAT V (m/s) 1.0E4 335.9 223.9
A (m/s2) 55.2 57.4 38.3
SIGQ (EKF) (m2 / s) 196.0 196.0 196.0
SIGQ (SDREF) (m2 /s4) 1.0E6 1.0E3 1.0E3
Meas Noise (%, im, tm) 0.1 50.0 50.0
Figures 17 and 18 show the positionstate response of the EKF and SDREF to a 1
percent initial error. (Note the oscillations in the SDREF estimate during the initial
settling portion of the flight.) Velocity and acceleration trends were similar.
The filters were then run with a 5percent initial error. (The initialization
parameters for the EKF and SDREF are given in Table 3.) Because of the type of
engagement and large distances and closing velocities involved, a fixedpercentage type
of initial error induces large disturbances in the xdirection. Therefore, little attention is
given to the performance in this direction, as these magnitudes of error are not likely to
occur. Figures 19 through 20 show the performance for a 5percent initial error. The
inability to properly adjust the initial covariance matrix in the SDREF causes a large
oscillatory error in the early portion of the estimation problem. This error was
unacceptable and a new initialization procedure was investigated.
The basic problem with the initialization of the SDREF is that the initial
covariance matrix is computed using the solution of the algebraic Ricatti equation. This
indicates that no large errors are expected. Therefore, there is no method of supplying the
SDREF with an initial covariance error of the magnitude of the expected initialization
errors. Because of this the SDREF computes a smaller magnitude covariance matrix and
therefore a smaller gain to be used in the state update. This is why the filter oscillates
during the initial portion of the flight. If the SDREF could be "jumpstarted" with a high
initial covariance, the majority of the initialization error could be taken out in the first
time step.
67
PPL
0 2 4 6 8 10 12
F1
0 2 4 8 10 1
Time (sec)
Figure 17. EKF position errors with 3ustandard deviations
for 1percent initialization error.
2.
1000
0 2 4 8 10
0 2 4 6 0 10
4o0
20
0 2 4 8
Time (sec)
Figure 18. SDREF position errors with 3astandard deviations
for 1percent initialization error.

~YU
INU
 111~1
II 10
E400
01s
1)0
so
S o
g 0
150
.60
40
S so
100
F
o 0 4 5 a 10
04 6 0
Time (sec)
Figure 20. SDREF position errors with 13standard deviations
for 5percent initialization error.
I I I
I
~
uD
E 'm
m
0 0 I 0 12
Io
o 4 e S 10
0 0 0 10 12
Time (sec)
Figure 19. EKF position errors with 3cstandard deviations
for 5percent initialization error.
^
Table 3. EKF and SDREF parameters fivepercent initialization error
Parameters (units) X Y Z
Aspect Angle (degrees) 30 30 30
SIGP (m) 5.0E3 168.0 112.0
SIGVP (m/s) 5.0E2 16.8 11.2
C2 (m2) 25.0E6 3.0E4 3.0E4
SIGXHAT o2 (m2 / 2) 25.0E4 3.0E2 3.0E2
CA2 (m2 / 4) 1.OEl 1.LEl 1.OEI
R (m) 1.0E5 10.0 10.0
XHAT V (m/s) 1.0E4 335.9 223.9
A (m/s2) 55.2 57.4 38.2
SIGQ (EKF) (m2/s4) 196.0 196.0 196.0
SIGQ (SDREF) (m2 / s4) 1.0E6 1.0E3 1.0E3
Meas Noise (%,am,Aum) 0.1 50.0 50.0
If the gain calculation used in both the EKF and SDREF is examined,
K= = p (_)Hk[H p, (_)H + R,]
it can be seen that the covariance used at time step zero, Pk (), is needed in the
calculation. If it is assumed that the value of Pk that the SDREF computes is Pk(+) (i.e.,
after the updated measurement has occurred), then the SDREF is initialized with P, () to
be used at time step zero. Using this methodology, the SDREF was initialized with the
same covariance as the EKF and the SDREF process noise was also reduced to the same
level as that of the EKF. This change in methodology had positive impacts on the
performance of the SDREF when confronted with initialization errors. Figures 21
through 26 show the new performance of the SDREF using the same initialization errors
as in the baseline runs. It can be seen that the initial oscillatory motion in the state
estimate errors has been corrected and the filter compares favorably with that of the EKF.
"I I i
Time (sec)
Figure 21. SDREF position errors with 3astandard deviations
after new initialization method, 1percent error.
I I
Time (sec)
Figure 22. SDREF velocity errors with 3oastandard deviations
after new initialization method, 1percent error.
Time (sec)
Figure 23. SDREF acceleration errors with 3astandard deviations
after new initialization method, 1percent error.
Time (sec)
Figure 24. SDREF position errors with 3cstandard deviations
after new initialization method, 5percent error.
I: I
It o I I
Time (sec)
Figure 25. SDREF velocity errors with 3astandard deviations
after new initialization method, 5percent error.
j ^I I
ITime (s )
 I ... m I Jc
Time (sec)
Figure 26. SDREF acceleration errors with 3ostandard deviations
after new initialization method, 5percent error.
i
This new initialization methodology not only improves the performance in the presence
of initial errors, but also allows a convenient way of initializing the SDREF.
6.3.2 Measurement Errors
This section investigates the performance of the EKF and SDREF when subjected
to various degrees of measurement noise. The relative performance of the two filters
using the baseline measurement noise was shown in Figures 3 through 5. Using the same
baseline initial conditions given in Table 1, two other levels of measurement noise were
added (Tables 4 and 5).
Table 4. EKF and SDREF measurement noise parameters level 1
Parameters X Y Z
Measurement Noise (%,rad,rad) 1.0 100.0E6 100.0E6
Filter Noise Variance (.01* Range)2 1.0E8 1.0E8
Table 5. EKF and SDREF measurement noise parameters level 2
Parameters X Y Z
Measurement Noise (%,rad, rad) 10.0 500.0E6 500.0E6
Filter Noise Variance (0.1 Range)2 2.5E7 2.5E7
The second level of noise is large for a typical engagement; however, it exercises
the robustness of each filter. Both filters were given measurement noise covariance
matrices that matched the noise characteristics given to the sensor. The position state
estimate plots are given in Figures 27 through 30. There is not a significant difference in
filter performance using the large noises. As expected, the estimates for both filters are
noisier and more erratic. The acceleration state estimate has become even worse under
these conditions. This estimate was initially poor, and with noisy measurements the
filters do not have enough information to estimate the target acceleration to the desired
degree of accuracy. It can also be noted that the SDREF estimation error under the
second level of noise seems to be worse than the EKF estimate. While these figures only
represent one run, the SDREF seems capable of adequate performance even under
stressing noise conditions.
6.3.3 Target Maneuvers
This section investigates the openloop performance of the EKF and SDREF
against two stressing target maneuvers. The capabilities of the filters against a thrust cut
off maneuver have been previously discussed. In this section two other nonlinear
maneuvers were tested. The first is a 10degrees per second rotating target maneuver.
This maneuver allows the target acceleration vector to rotate approximately 110 degrees
during the engagement. This is a moderately high target maneuver, but it is possible for
an ICBM to perform this type of action in an active countermeasure scenario. The second
stressing engagement consists of a dogleg maneuver in which the target acceleration
vector rotates 90 degrees out of plane and then reverses itself to rotate back inward 90
2OO
300
soo
10
x 100
BOO
30
20
3O
10
i o
N 0
30
Time (sec)
Figure 27. EKF position errors with 3oastandard deviations
1percent error in X, 100 microradian error in Y and Z.
Time (sec)
Figure 28. SDREF position errors with 3ostandard deviations
1percent error in X, 100 microradian error in Y and Z.
I ?A..S
I I
t0 0 iiiir i 
sQ    .
mo=^^
E m ;  
a 2  *l .
S. .ao ,' 
x ) : == = 
o  
.80  
.100 *_ '
. ............ ..
0 2 10 1
so
2w
015
,0 2 ^     4 S 1
Time (sec)
Figure 29. EKF position errors with +3" standard deviations
10percent error in X, 500 microradian error in Y and Z.
0 2 4 6 6 10 1o
Time (sec)
Figure 30. SDREF position errors with 3astandard deviations
10percent error in X, 500 microradian error in Y and Z.
E
1 IcT I 1
a ra 12
76
degrees. The initial conditions used for this analysis were the same as the baseline
scenario except for the process noise characteristics given to the filters (see Table 6).
Table 6. Maneuvering target parameters
Parameters (units) X Y Z
Aspect Angle (degrees) 30 30 30
SIGP (m) 0 0 0
SIGVP (m/s) 0 0 0
o2 (m2) 1.0E5 1.0E5 1.0E5
SIGXHAT o2 (m2 / 2) 1.0E3 1.0E3 1.0E3
Cr2 (m2 / s4) 1.0E2 1.0E2 1.0E2
R (m) 1.0E5 10.0 10.0
XHAT V (m / s) 1.0E4 335.9 223.9
A (m/s2) 55.2 57.4 38.2
SIGQ (EKF) (m2 /s4) 866.0 866.0 866.0
SIGQ (SDREF) ( / s4) 1.0E3 1.0E3 1E3.0
Meas Noise (%,nm, rn) 0.1 50.0 50.0
The process noise levels were tuned to allow the filters to open up and "accept" the
unexpected target maneuvers. While this allows for better tracking, it does lead to a
noisier state estimate in both filters. Figures 31 through 34 show the relative performance
of each filter against these two scenarios.
As can be seen from these figures, there are no appreciable differences in the state
estimates between the two filters. They each track this maneuverable target, however the
state estimate error is noisier and larger, as expected. It should be noted that once the
initial transient terms in the EKF covariance have settled, they are similar in shape and
value to the "steadystate" covariance terms of the SDREF.
45 ..8 0 12
0 2 4 6 S 10 02
Time (sec)
Figure 31. EKF position errors with 3astandard deviations
for 10 deg/sec rotating target maneuver.
 20
0 2 4 0 6 10
... .. _
0 4 6 
Time (sec)
Figure 32. SDREF position errors with 3ostandard deviations
for 10 deg/sec rotating target maneuver.
.10 4 10 12____________
B8 10
lo
lo
0 0 4 4 8 10 12
o 04 6 l 0
Time (sec)
Figure 33. EKF position errors with 30 standard deviations
for dogleg target maneuver.
8 o
. Ii K ;. / , Avv '. A
.1
Time (sec)
Figure 34. SDREF position errors with 3+o standard deviations
for dogleg target maneuver.
IL1
w~
~nu
iraaoa
I to
79
6.4 SDREF Closed Loop Performance
6.4.1 Other Engagement Aspect Angles
Up to this point the EKF and SDREF performance in an openloop environment
have been examined. The following analysis examines the closedloop performance of
each filter against various engagement scenarios. This includes a variety of engagement
aspect angles, target maneuvers, initial condition errors, and measurement noise levels.
In this section, however, performance in terms of probability of hit and percentage of fuel
used to achieve the hit are measured. Measuring the percentage of fuel used, allows a
comparison of how noisy the state estimates are, as increased error in the beginning of the
engagement results in fuel being wasted performing unnecessary maneuvers. This
analysis was performed using 31 Monte Carlo run statistics. The first analysis focused on
three engagement scenarios. Up to this point a 30,30,30degree closing aspect angle was
used. A headon, 90degree beam shot, and a tailchase engagement are now added.
Figure 35 shows the three engagements and the corresponding simulation parameters are
given in Tables 7 through 9.
Y
(0,0,0) Aspect Angle Head On
(0,0,0) Aspect Angle Head On
13*
(0,0,90) Aspect Angle Beam Shot
Y
(0,150,0) Aspect Angle Tail Chase
Figure 35. Aspect angle description.
x r
t
Table 7. Headon simulation parameters
Parameters (units) X Y Z
Aspect Angle (degrees) 30 30 30
SIGP (m) 0 0 0
SIGVP (m/s) 0 0 0
r2 (m2) 1.0E5 1.0E5 1.0E5
SIGXHAT O2 (m2 / s2) 1.0E3 1.0E3 1.0E3
oaA (m2/s4) 1.0E2 1.0E2 1.0E2
R (m) 1.0E5 10.0 10.0
XHAT V (m/s) 10000.0 335.9 223.9
A (m/s2) 55.2 57.4 38.2
SIGQ (EKF) (m2 / s4) 866.0 866.0 866.0
SIGQ (SDREF) (m2 /s4) 1.0E3 1.0E3 1.0E3
Meas Noise (%,/pm, mn) 0.1 50.0 50.0
Table 8. 90degree beam shot simulation parameters
Parameters (units) X Y Z
Aspect Angle (degrees) 0 0 90
SIGP (m) 0 0 0
SIGVP (m/s) 0 0 0
62 (m2) 1.0E5 1.0E5 1.0E5
SIGXHAT O,2 (m /s2) 1.0E3 1.0E3 1.0E3
___ (m /s4) 1.0E2 1.0E2 1.E2
R (m) 1.0E5 10.0 10.0
XHAT V (m/s) 1.0E4 498.2 0.0
A (m/s2) 0.0 88.3 0.0
SIGQ (EKF) (m2/s4) 196.0 196.0 196.0
SIGQ (SDREF) (m2 /4) 196.0 196.0 196.0
Meas Noise (%,mn, Um) 0.1 50.0 50.0
82
Table 9. Tailchase simulation parameters
Parameters (units) X Y Z
Aspect Angle (degrees) 0 150 0
SIGP (m) 0 0 0
SIGVP (m/s) 0 0 0
a,2 (m2) 1.0E5 1.0E5 1.0E5
SIGXHAT O,2 (m / s2) 1.0E3 1.0E3 1.0E3
o2 (m2 / s) 1.0E2 1.0E2 1.0E2
R (m) 1.0E5 10.0 10.0
XHAT V (m/s) 1.0E4 0.0 237.9
A (m/s2) 76.5 0.0 44.1
SIGQ (EKF) (m2 /s4) 196.0 196.0 196.0
SIGQ (SDREF) (mi / s4) 196.0 196.0 196.0
Meas Noise (%,/Um, pm) 0.1 50.0 50.0
The 90degree beam shot is the most stressing of these engagements as the entire
target acceleration profile is out of plane. The headon case is one of the easier
engagements when using a range measurement, and the tail chase lies inbetween these
two in terms of estimation difficulty. Table 10 displays the results from the Monte Carlo
analysis against these engagements. Probability of hit (PHIT) is the probability that the
interceptor hit the target. Circular error probable (CEP) is the probability that 50percent
of the misses lie within this bound, 95percent miss includes ninetyfive percent of the
misses, and worst miss displays the worst of the Monte Carlo runs. The mass usage
numbers reflect the same quantities except in terms of percent of divert fuel used. The
vehicle has a mass of 4 kilograms, 1.5 of which is devoted to divert fuel. When this limit
is exceeded, the interceptor will continue to fly towards the target with no further
guidance corrections, typically resulting in large miss distances. Since these engagements
had good initialstate estimates and the target is not maneuvering, all of the runs were
hits; however as shown in the following sections that this is not always the case.
Table 10. Miss distance performance for varying scenarios
Parameter (units) 00 90 1500
EKF SDREF EKF SDREF EKF SDREF
PHIT 1.0 1.0 1.0 1.0 1.0 1.0
CEP (m) 0.185 0.262 0.432 0.625 0.315 0.265
Mean Miss (m) 0.213 0.272 0.546 0.667 0.347 0.354
95% Miss (m) 0.364 0.426 1.212 1.209 0.654 0.671
Worst Miss (m) 0.421 0.599 1.828 1.679 0.803 0.895
Mass CEP (%) 16.70 27.61 67.38 72.67 43.42 48.06
Mass Mean (%) 16.90 28.27 66.98 72.76 43.11 48.40
Mass 95% (%) 21.03 32.78 69.25 75.22 46.10 51.62
Mass Worst (%) 22.70 33.16 71.20 75.51 48.57 52.59
From Table 10, no significant difference in the closedloop performance of these two
filters under the above scenarios is shown. All performed adequately and ensured a hit
against the target in all cases. Fuel usage was also not a decisive factor as both filters
differed by only a few percentage points.
6.4.2 Stressing Target Maneuvers
To further test filter performance, two other target maneuvers were considered.
These maneuvers are on the outer edge of filter capability of properly tracking the "true"
state. They were included to investigate whether the parameterization of the SDREF
allows an advantage over the linearization of the EKF in highly nonlinear environments.
The first maneuver is a rotating target, however in this engagement the engagement
aspect angle was a 90degree beam shot and the rotation rate was increased to 20 degrees
per second. Therefore, the target will have completed an approximately 200degree turn
during the length of the engagement.
The second maneuver consisted of a dogleg evasive maneuver. This maneuver is
similar to the previously described dogleg maneuver except that it was performed in the
90degreebeamshot scenario, causing it to become more difficult for the filter to
estimate.
Table 11 displays the results of Monte Carlo runs. The process noise levels for
each filter were tuned to give the best performance against these scenarios.
Table 11. Miss distance performance for stressing target maneuvers
Parameters (units) 200/Sec Dogleg
EKF SDREF EKF SDREF
PHIT 0.935 0.806 0.871 0.742
CEP (m) 1.61 2.20 1.14 2.85
Mean Miss (m) 1.83 200.20 1.50 3.05
95% Miss (m) 4.30 7.44 3.94 6.95
Worst Miss (m) 6.00 6121.89 5.59 7.36
Mass CEP (%) 100 100 100 100
Mass Mean (%) 99.94 99.99 99.87 99.98
Mass 95% (%) 100 100 100 100
Mass Worst (%) 100 100 100 100
As shown in Table 11, the SDREF performed slightly worse against the targets than the
EKF. The one large miss distance for the SDREF was caused by poor tracking, which
caused the interceptor to run out of fuel early and coast to a large miss distance.
The purpose of this section was to stress the filters to their extreme operational
limits and investigate the behavior. These types of target maneuvers would be considered
85
unlikely. However, they do provide a convenient way to determine the robustness of each
filter design. Even at these extremes, the SDREF was almost able to meet the EKF
performance measures.
6.4.3 Measurement Errors
This paragraph discusses filter sensitivity to various levels of measurement error.
For this analysis the range measurement error was held constant at 1 percent due to the
size of the range error when dealing with large engagement distances. Azimuth and
elevation noises were included for 100, 500, and 1000 microradian errors. The results
displayed in Table 12 were done using a 31run Monte Carlo experiment.
Table 12. Miss distance performance for varying measurement noise
Parameters (units) 00rad 500rad 1000rad
EKF SDREF EKF SDREF EKF SDREF
PHIT 1.0 1.0 0.968 1.0 0.935 0.806
CEP (m) 0.455 0.460 1.39 1.41 3.09 2.22
Mean Miss (m) 0.468 0.568 1.47 1.45 3.17 3.14
95% Miss (m) 0.775 1.218 2.47 2.60 5.66 6.29
Worst Miss (m) 0.935 1.569 2.55 3.08 6.03 8.18
Mass CEP (%) 77.13 76.86 84.19 81.83 89.79 86.99
Mass Mean (%) 77.57 76.96 84.66 82.36 90.23 87.05
Mass 95% (%) 80.67 78.95 88.27 85.12 94.46 89.40
Mass Worst (%) 81.23 79.55 90.92 86.31 96.66 92.67
Once again, the results between the two filters are similar. Even with highmeasurement
noise levels, both filters performed adequately, missing only one or two times in the
Monte Carlo run. The fuel usage numbers are also comparable. These results verify the
results of the pdf comparison between the filters, as the probability distributions are
similar between the EKF and SDREF. Therefore, it can be assumed that the filters are
Monte Carlo run. The fuel usage numbers are also comparable. These results verify the
results of the pdf comparison between the filters, as the probability distributions are
similar between the EKF and SDREF. Therefore, it can be assumed that the filters are
using approximately the same set of distribution functions, even though they are achieved
in different ways.
6.4.4 Initialization Errors
Using the new initialization method described in section 6.3.1, a closedloop
analysis of the EKF and SDREF was performed. The initial conditions used for the open
loop initialization study were used for this investigation. The filters were run closedloop
in a 31run Monte Carlo simulation. Table 13 displays the results of these runs. Using
this new initialization method not only shows better performance in the singlerun, open
loop case, but performs equally as well in the closedloop analysis. Probability of hit,
miss distance, and fuelusage numbers show that the SDREF performance is similar to
that of the EKF.
Table 13 Miss distance performance for varying initialization error
1% Error 5% Error
Parameters (units) E SDREF EKF5% Errr
PHIT 1.0 1.0 1.0 1.0
CEP (m) 0.436 0.377 0.336 0.389
Mean Miss (m) 0.464 0.380 0.423 0.432
95% Miss (m) 0.849 0.750 0.858 0.770
Worst Miss (m) 1.01 0.793 1.13 1.12
Mass CEP (%) 75.13 79.18 72.53 92.03
Mass Mean (%) 75.00 79.45 72.53 92.03
Mass 95% (%) 77.34 81.92 74.56 93.20
Mass Worst (%) 78.12 84.78 75.55 93.42
87
6.5 Summary of Exoatmospheric Guidance Problem
The exoatmospheric guidance problem analysis demonstrated, with some minor
caveats, the SDREF was capable of satisfying the interceptor guidance problem. The
assumptions of detectability in solving the SDRE required the addition of a range
measurement. An ad hoc initialization algorithm needed to be developed to properly
initialize the SDREF. With these two restrictions in place, however, the SDREF
performed similarly to the EKF for the exoatmospheric guidance problem for several
engagements, noise sources, and initialerror combinations.
CHAPTER 7
PENDULUM PROBLEM FILTER DEVELOPMENT
To evaluate the effect of nonlinear dynamics on the SDREF performance, a simple
twodimensional pendulum problem was constructed. The scenario consists of a simple
rigid pendulum with friction about the connection point, see (Figure 36). The
measurements are assumed to come from an accelerometer attached to the bob, measuring
g.
\0
/.
Figure 36. Pendulum setup.
89
The pendulum is given some initial position and velocity and the friction term
damps the oscillations. This system, while simple, demonstrates nonlinear dynamics and
nonlinear measurements together in the same problem.
The system dynamics for the pendulum problem are given by the 2"d order
nonlinear differential equation
B=gsinO BO (7.1)
L
where BO is the friction term. B was set to be 0.5.
Both the SDREF and EKF were derived using the stochastic nonlinear system
.i = f(x) +w
z= h(x)+v
where w and v are Gaussian, zeromean whiteprocess and measurement noise,
respectively, given by
E[w(t)w'(t + r)= W(t)3(r)]
(7.3)
E[v(t)v(t + r)= R(t)S()]
If
x, =
x (7.4)
then the state equations for the pendulum system are
xi = x
g (7.5)
2 = in(x,)Bx2 (7
L
The measurement equation is
z = sin(x1) Bx2
L

Full Text 
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID E3A4DIPYE_PKLNE8 INGEST_TIME 20130123T15:27:35Z PACKAGE AA00012992_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES
PAGE 1
AN ANALYSIS OF A NEW NONLINEAR ESTIMATION TECHNIQUE: THE ST ATEDEPENDENT RICATTI EQUATION METHOD By CRAIG M. EWING A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL ......u.. .L.L..il.J. MENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1999 ..
PAGE 2
ACKNOWLEDGMENTS I would like to thank my advisor, Dr. Norman FitzCoy, for hi s efforts in co ntinually pushing me to develop a better dissertation. I agree it is now someth in g I can look upon with pride. I also thank Dr. Jame s Cloutier, the man behind the idea Withou t his help and guidance this work could never have been written. I would especially like to thank my wife, Darsi, for her neverending s upport through these years of work. She never let me give up 11
PAGE 3
TABLE OF CONTENTS page ACKNOWLEDG"MENTS . ... ... . .... .... .... .... I ii LIST OF TAB LES .... ..................................... .................. .... ....... ......... .................. .... .... v LIST OF FI GuRES ... .... ....... . ..... ............ .... .......... ....... ....... .... ........ .... .... .... .................. vi AB ST .RA CT ... ..... ... .... ....... ..... . . .... .... . . .... .......... ...... .... ....... ................ ............................. X CHAPTERS 1 'IN"TRODUCTION ............. . .................... ....... ..... ..... .............................. .... ............ .. 1 2 LITERATURE SUR "VEY ................. .... .... ............. .... ...... .... .... .... .... . ...................... .... 4 3 SCOPE OF WORK ................... ........... . ......... ....... ... ... ............. .............................................. 8 4 EXOATMOSPHERICGUIDANCEPROBLEM ENGAGEMENT SCENARIO ..... 10 5 EXOATMOSPHERICGUIDANCE PROBLEM FILTER DE"VELOP"MENT ... ...... 13 5. 1 EKF Deri vation .. . . . . .. . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. .. . . .. .. . .. . . . .. .. .. .. . .. .. .. .. .. . . . 14 5 .2 Boot s trap Estimator Derivation ....... ....... ....... ............... .... ....... . ....... . .................. 19 5.3 SDREF Derivation ............. . ....... .... ......................... . ...... ....... ........................ ........... ...... 22 6 EXOA TMOSPHERIC GUIDANCEP ROBLEM SIMULATION RESULTS .... ....... 49 6.1 Baseline Performance ................................. ................. ......... .................... .... ...... . ....... ........... 49 6.2 PDF Distribution Comparison ............. ................ ....... .... ........................ ......... .... 54 6 .. 3 Sensitivity Analysi s ........ .............................................................. ........... ........... ................................. 64 6.4 SD REF Closed Loop Perforrnanc e ................................ . .... .......................... ....... 79 6 5 Summary of Exoatmospheric Guidance P roblem .... ......................... .......... ..... ..... 87 7 PENDUL UM PROBLEM FILTER DE"VELOPMENT .... ............ . ....... ............... ....... 88 7. I SD REF Derivation ......... .... .......... .............................................. .......... ................... ... ... 90 7 .2 EKF Development .... ...................... ....... ...... ...... ...... ........ ....................... .................... 91 8 PENDULUM PROBLEM SIMULATION RES UL TS . .... ............ . ..... . ....... ...... ..... 92 111
PAGE 4
9 MSDREF DEVELOPMENT ...... .... ......... ..... ...... ... ......................... . .................. .... 103 9 .1 MSD REF Derivation ... . .......... .... ... ......... .... .... ..... . .... ............ .... .... ..... ......... 103 9 .2 Stability Proof ......... ............... ....... ... ........ .... ............. . .... ......... .... .... ................. 104 10 MSDREF PROBLEM SIMULATION RES UL TS ......................... ......... ... .......... 112 11 CRITICAL ANALYSIS OF RESULTS ....... ....... ................ ... ............ ................. ... 119 REFERENCES ................................................... .. ............ ....... . ... ... ........ ....... ......... .... ........................ 122 BIOGRA.PHICA L SKETCH ........................................ .. .......................................... ........ ............................ ........... .... ............. 125 IV
PAGE 5
LIST OFT ABLES Table page 1. Filter simulation parameters baseline conditions ..... ... . .... ...... .... .... . ........ .... ... ...... ... 50 2. EKF and SD REF parameters onepercent initialization error . .................................. 65 3. EKF and SDREF parameters fivepercent initialization error ....................... ..... ...... 69 4. EKF and SD REF measurement noise parameter s level 1 ................. ..................... 72 5 EKF and SD REF measurement noise parameters l evel 2 ..................................... . 72 6. Maneuvering target parameters ......... . ..... . .... ......... ........................ ....... .... ............ 76 7. Head on simulation parameters .................. ..... ......... .... . ..... ... ............... ...... . ........ .... 81 8. 90degree beam shot simulation parameters ............................... .... .... .... ................... 81 9. Tailcha se sim11lation parameters ....... ..... . .............................. ......... ...................... 82 10. Miss distance perforrnance for varying sce narios ... . ....... .... .... ..... . .................. ...... 83 11. Miss distance perforrr1ance for stressing target maneuvers ............ .... .... .... ........... 84 12. Miss distance performance for varying measurement noise . ............................... .... 85 13. Miss distance performance for varying initialization error .............................. .... ...... 86 V
PAGE 6
LIST OF FIGURES Figur e p age 1. Generic intercept vehic l e .... . ..... .... .... . ......................... ............... ......... .... .... ......... 11 2. Guidance system block diagram .... .... ...... ....... . .... ..... ... ............ ..... . ... .... . ........ ..... 12 3. EKF pos ition error with ostandard deviations ... ............ . .... ............ ....... . . ........ 52 4. SD REF position error with ostandard deviations ........ . ........ . ....................... .... 52 5. Bootstrap posi tion error with os tandard deviations ....... ............... ..... . .... ............ 52 6. Comparison of EKF and SDREF covaria nce os tandard deviation s for Y po s ition ........ . .... .................... ............... ........ .......... .... ............... ....... .......... 53 7. Expanded view of Figure 6 ..... ... .......... ....... .......... ..................... ........................ . 53 8 Po ste rior probability density functions at time= 1 5 and 9 seconds for X position ............... .... . ..... ............... .... .... ..... ........... .... . .................... ......... 55 9. Pos terior probability density function at time= 1 5, and 9 seconds for Y position .... . .... ... . ......... .......... ....... . ...... ... ... .... ................ . ...... ... ... ... ...... ..... 56 10. Monte Carlo pos terior probability density functions at time= 1 5 and 9 seconds for X position ....... . ..... .... ....... ..... .... ........ . .................. .... . ............. .... 58 11. Monte Carlo poste rior prob a bility d e n s ity functions at time= 1 5 and 9 seconds for Y position ....................... ................... .......... ......... .... .... ............ . .... 59 12 EKF estimated and true acceleration . ... ........ .... ............. ..... .... .... ........... ....... ......... 61 13. SD REF estimated and true acceleration ............... ....... ..... ....... .... ........................... 61 14 Bootstrap estimated and true acceleration ... ................................ .... . ................. ...... 61 15. Monte Carlo posterior probability den s ity functions at time= 2, 2.5, and 3 seconds for thru s t cutoff X position ........ . ......... ............. . . .... ................... . ..... . 62 VI
PAGE 7
16. Monte Carlo posterior probability density functions at time= 2, 2.5, and 3 seconds for thrust cutoff Y position ...... ..... ................... . ....... .... .... ........... ..... 63 17. EKF position errors with a standard deviations for I percent initialization errors ..... ............ ...... ........................................... .................... ........ 67 18. SDREF position errors with standard deviations for Ipercent initialization enors .................. .... .................. ... ... .................. ............... . ........... 6 7 19. EKF position errors with standard deviations for 5percent initialization errors ........ .... .................. ...................... ...... ................. .................... 68 20. SDREF position errors with a standard deviations for 5percent initialization errors ... ....... ............... ... ... ............... . .... .................. . ............ .... ..... 68 21. SDREF position errors with a standard deviations after new initialization method, Ipercent error ... . ......... .... ...................... .... ..................... . ............ ..... .... 70 22. SDREF velocity errors with a standard deviations after new initialization method, 1percent error .... ............... ....... ............... .... .......... ..... . ......... . .... . .......... 70 23 SDREF acceleration errors with a standard deviations after new initialization method, 1percent error .............................. ........ ............ ..... ........... .... .................. 7 0 24. SDREF position errors with a standard deviations after new initialization method, 5percent error .... ....... . ............. .... ...... ........ ............... .... ......... ..... ......... 71 25. SDREF velocity errors with a standard deviations after new initialization method, 5percent error ....... ....... ........ . .... ................ . .......... ............ ................. 71 26. SDREF acceleration errors with a standard deviations after new initialization method, 5percent error ........................................... . .... ......................................... 71 27. EKF position errors with a standard deviations Ipercent error in X, 100 microradian error in Y and Z .......... ................ ............................... ............... 7 4 28. SDREF position erro rs with a sta ndard deviations Ipercent error in X, 100 microradian error in Y and Z ...... ....... ............. ............................... ............. ... 7 4 29. EKF position errors with a standard deviations 10percent error in X, 500 microradian error in Y and Z ...................... .......... ............. ...... .... .......... ....... 75 30. SDREF position errors with a standard deviations IOpercent error in X, 500 microradian error in Y and Z ....... .... ........ . .... ............. .... .... .... ..................... 7 5 Vll
PAGE 8
31. EKF position errors with a sta ndard deviation s for 10 de g/sec rotating targe. t maneuver .... ......... .... ....... .... .... ..................................... .... . ....... .... ....... ..... 77 32. SDREF pos ition errors with a sta ndard de vi ation s for 10 d eg/sec rotating target maneuver ........................................................................... ......... .......... ...... 77 33. EKF pos ition errors with a s tandard deviations for dogleg target man euver ...... 78 34. SDREF pos ition errors with a standard deviation s for dogleg target maneu ve r 78 3 5 A s pect angle de scr iption ............................................ .... .... .... .... . . . ......... . ........... 80 36. Pendulum setup . ... ... .... .... .... ...................................... .... ...................... . ......... ...... 88 37. EKF and SD REF estimates vs. true theta with no initiali zatio n e rror ... .... . . . ... ..... 93 38. EKF and SD REF est imate s vs. true thetado t with no initialization error ................ 93 39. EKF and SDREF estimation error for theta with no initialization error ... .... ... ......... 94 40. EKF and SD REF esti mation error for th etado t with no initialization error .... .... .... 94 41. Monte Carlo probability d e n s ity function for theta with no initialization error, time= 0.5 sec ... ....... ..... ................... .. .. .... .............. ............ ............... ... ..... ........ .... 95 42. Monte Carlo probability den si ty function for theta dot with no initiali z ation error, time = 0. 5 sec ............ .... ........... ............ ..... ....... ..... ........ .... ....................... 9 5 43. Monte Carlo probability d e n s ity fun ctio n for theta, with no initialization erro r time= 2.0 sec ... ................. .......... .... ..... ................... ........... .... ........................ 96 44. Monte Carlo pr oba bility d e n s ity function for theta dot with no initialization error, time = 2 0 sec ...... ................. ........... . ..... .... ................... . .... .... .................. 96 45. Monte Carlo probability density function for theta with no initialization e rror, time = 3. 0 sec ... . .... ............ ..... .... ............ .... ........... ..... .... .... ............... ............ 97 46 Mont e Carlo probability d ensi ty function for thetadot, with no initialization error, time = 3 .0 sec .... ....... .... ....... . ....... .... .......... .......... ... ..... .... . . ............... ... 97 47. EKF and SD REF theta estimation error with initialization error . . .... ....... .............. 99 48 EKF and SD REF thetadot estimation error with initialization error .................... ... 99 49. Mont e Carlo probability density function for theta, with initiali zatio n error, time = 0 5 sec . . . ......................... .... . . ..... .... ................ .... .... .... ............ ......... 100 Vlll
PAGE 9
50 Monte Carlo prob a bility density function for thetadot, with initialization eno1, time= 0.5 sec ................ .......................... ...... ...... ............................ ............................ .. .......... .. ...................... .. ..................... ..... 100 51. Monte Carlo probability density function for theta with initialization error, time= 2.0 sec . .... . ....... .................... . ... .. ...... .... ... ... ... ... .... . ................... ......... 101 52. Monte Carlo probability den s ity function for thet adot with initialization e rror tim e = 2.0 sec ................................................... ........ .................................................................................................. ................... 101 53. Monte Carlo probability d e n s ity f unction for theta with initializ a tion e rr o r time= 3.0 sec ..................................................................................................... ................................................................ .. ........................ 102 54. Monte Carlo probability den s ity function for thetadot, with initiali za tion e rr o r time = 3 0 sec .................... .... .... .... .... .... ............. .... ................. ... ......... . .......... . 102 55. True vs. estimated Xl state ... .......... .... .... ........................... .... .... .... ...... ....... . ....... 117 56. True vs estimated X2 state ........... .... ....... ...... .... . .... .... ........ ..... . ... . .............. .... 117 57. True vs. estimated X 1 state, nonzer o initial conditions . .... . .... . . .... . ....... . . . .... 118 58 True vs. estimated X2 s tate non ze ro initial conditions ... .... ....... . . . . .... ..... ......... 118 IX
PAGE 10
Abstract of Dissertation Presented to the Graduate School of The University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy AN ANALYSIS OF A NEW NONLINEAR ESTilvlATION TECHNIQUE: THE STATEDEPENDENT RICATTI EQUATION METHOD Chairman: Norman FitzCoy By Craig M. Ewing August 1999 Major Department: Aerospace Engineering, Mechanics, and Engineering Science Research into nonlinear estimation techniques for terminal homing missiles has been conducted for many decades. The terminal state estimator, also called the guidance filter, is responsible for providing accurate estimates of target motion for use in guiding the missile to a collision course with the target. Some form of the extendedKalman filter (EKF) has become the standard estimation technique employed in most mode111 weapon guidance systems. EKF lin earization of nonlinear dynamics and/or measurements can cause problems of divergence when confronted by highly nonlinear conditions. The objective of this dissertation is to analyze a new nonlinear estimation technique that is based on the parameterization of the nonlinearities. This parameterization converts the nonlinear estimation problem into the form of a steadystate continuous Kalman filtering problem with statedependent coefficients X
PAGE 11
This new technique, called the s tatedependent Ricatti equation filter (S DREF), allows the nonlinearities of the system to be fully incorporated into the filter design, before stochastic uncertainties are imposed, without the need for linearization. The SDREF was investigated in three problems : an exoatmospheric, terminal homing, ballistic miss ile intercept problem; a highly nonlinear pendulum example; and an algorithmic lo ss of observability problem. The exoatmospheric guidance problem examined nonlinear measurement s with linear dynamics To investigate the SDREF when used with a combination of nonlinear dynamic s and nonlinear measurement s, a highly nonlinear, twostate pendulum problem was also examined. While the se problems were useful in gaining insight into the performance characteristics of the SDREF, no formal proof of stability could be determined for th e original formulation of the estimator. The original SDREF so lved an algebraic SDRE that arose from an infinite time hori zo n formulation of the nonlinear filtering problem. A modification to the SDREF formulation was developed that led to a differ e ntial SDREF and a proof for local asymptotic stability was achieved. The modification r e moved the infinite time horizon assumption and integrated the coupled s tatedependent state and covariance equations. Thi s new form of the e s timator is called the modified SDREF ( MSDREF ) A problem involving algorithmic los s of observability was then examined. This problem shows a performance advantage when u s ing parameteri zat ion versus l i nearization as in the EKF technique. X l
PAGE 12
CHAPTER 1 INTRODUCTION The extended Kalman filter (EKF) bas become an accepted standard for both tactical and strategic interceptor guidance problems. While some problems lend themselve s to this type of estimation s trategy systems with highly nonlinear dynamic environments or misunderstood uncertainties can cause divergence of the filters. Many alternative nonlinear estimation techniques have been researched in an attempt to overcome this problem. The linearization of the nonlinear dynamics and/or measurements is one area of concern The EKF linearizes about the current estimate of t he state. If the state estimate is fairly accurate, thi s method has been shown to be effective. If, how ever, the estimates are poor, linearization can cause errors in the filter, which can lead to degraded performance or po ss ible divergence of the state estimate. It is this linearization that motivates the following research. This dissertation investigates a new nonlinear estimation algorithm, which originates from a recent nonlinear regulator design technique It uses parameterization to bring the nonlinear sys tem to a linearlike structure with statedependent coefficients ( SDC ). The filtering equivalent to the regulator technique can be realized in a recursive nonlinear estimation algorithm, which has the structure of the steadystate continuous Kalman filter. The Kalman gain is found by solving the sta tedependent Ricatti equation at each filter update time. The filtering technique, called the statedepende nt Ricatti equation filter (SDREF), I
PAGE 13
2 allows the incorporation of nonlinear dynamics and/or measurements in the filter design without the need to linearize the equations. The SDREF was analyzed in a several ways to determine the benefits and issues associated with this type of structure. By using a sampled data Bayesian filter, known a s the bootstrap filter the ''true'' probability density function (pdf) of the SDREF states was monitored. This allowed an investigation into the differences in pdf evolution for linearizing versus parameterizing filters. Finding the true pdf can be difficult; however recent work, by the author on the bootstrap estimator, proved valuable for this analysis. The bootstrap technique uses thousands of samples of the pdf along with Bayes' rule to estimate the true pdf of any given state With this filter and a large number of samples t he ''true'' state pdf can be closely approximated. An EKF was evaluated in a similar manner for comparison purposes The SDREF was exercised in three problems. A realistic and detailed evaluation of the filter was accomplished by examining an exoatmospheric guidance problem. The SDREF, bootstrap filter, and EKF were coded into a sixdegreeoffreedom (6DOF ) FORTRAN simulation for analysis in an exoatmospheric terminal homing engagement using passive and active sensor infortnation. The sensor measurements included azimuth elevation and range. The characteristics of the SDREF and EKF were evalt1ated in several different engagement scenarios. The SDREF was developed in a Cartesian coordinate frame, which causes the measurements to be nonlinear, while the dynamic s remain linear The second problem was that of a highly nonlinear, twodimensional pendulum problem Evaluation of the SDREF when used with nonlinear dynamics and nonlinear
PAGE 14
3 measurements was examined. The SDREF and EKF were coded into the visual simulation environment VisSim. The performance of the filters and their pdf' s were compared. This academic problem was well suited to exploring any differences that nonlinear dynamics may have on the performance or pdf characteristics of the SDREF. No forrnal proof of stability could be determined using the above formulations of the SDREF, so an alternate filter formulation was sought. The original SDREF was based on an infinitetime horizon formulation of the nonlinear filtering problem that resulted in an algebraic SDRE. A MSDREF was developed that utilizes a differential SDRE. Because the Ricatti equation is state dependent, it must be integrated in conjunction with the state estimates. A proof of local asymptotic stability was able to be determined. The new MSDREF was examined in a problem consisting of loss of algorithmic observability; and the results indicated a significant perfortnance advantage over the EKF. This dissertation fully develops the requirements for the SDREF and MSDREF. It includes an investigation of controllability and observability, which are pointwise and in the linear sense and pertain to a valid SDRE solution. An exhaustive investigation of the pdf characteristics, noise sensitivity, initialerror robustness, and closedloop performance of the SDREF using the EKF as a comparative source was performed A theoretical basis for the MSDREF, when used as an observer, was shown. The following research provides theoretical and practical infor1nation demonstrating that the SDREF and MSDREF are viable candidates for nonlinear estimation, and for some problems can out perform the EKF.
PAGE 15
CHAPTER2 LITERATURE SURVEY The literature abounds with research of various nonlinear estimation technique s, most of it mathematical and not application oriented. An elegant presentation of the hi st ory of airtoair target s tate estimation techniques by Cloutier et al.1 outlines a wide variety of techniques which include the Kalman filter, EKF, secondorder Gaussian fil ter (SOGF), recursive nonlinear filter recursive maximumlikelihood filter, and the assumed density filter. These filters have also been examined in the context of singleand multiplemodel adaptive filters. Adaptive filters may make use of certain propertie s of the residual or covariance to continuously change filter parameters to gain better filter performance On the other hand, multiplemodel filters use banks of filters to model many possible target dynamic situations. Various methods are then used to either select the results of one filter or combir1e that of seve1al. With the advent of the Strategic Defense Initiative ( SDI), a new set of problems arose when estimating the motion of a boosted intercontinental ballistic missile (ICBM) with the use of passiveonly measurements. While the EKF i s still the mos t widely accepted estimator for this type of engagement, it has been shown to have degraded performance when stressed by highly maneuverable targets. 2 4 This is the same problem that tactical filters experience in maintaining track of highly maneuverable aircraft or missiles. When tracking a highly nonlinear maneuver, state estimates can become poor and the linearization assumption begins to break down. The specific problems associated 4
PAGE 16
5 with guidance filters for the exoatrnospheric intercept problem have been addressed in several research programs noted in the following paragraph. Schmidt2 investigated an extension of work by Daum in which a closedforrr1 solution to the FokkerPlanck equations for propagating the pdf was developed A target maneuverdetection algorithm was added to this formulation to enhance performance The emphasis in this research was to more accurately model the target acceleration based 011 the rocket equation, instead of a Markovprocess. The target model caused the state equations to become nonlinear as well. The results of the work did show an improvement over the EKF technique; however this was accomplished, in part through the u se of the maneuverdetection algorithm Gido34 developed a nonlinear estimation algorithm, which numerically implements the theoretical work of Kolmogorov. The algorithm places a grid over the pdf and use s the FokkerPlanck equations coupled with Bayes rule to provide the state estimates. While this method showed improved performance over the EKF, it is computationally intensive. In the author's masters' thesis at the University of Florida ,5 another novel nonlinear estimation technique wa s investigated It involved using thousands of samples of the pdf, along with the state dynamics and Bayes' rule to propagate and update the pdf .6 7 The technique places no restrictions on the type of statistics used to model noise It al so makes no linearization assumptions. While computationally intensive, it is directly applicable to a massively parallel proces s ing architecture since the same operation is done to each sample. While the bootstrap method may not be seen in any missile guidance system design in the near future, it does provide a brute fo rce method for obtaining a reasonable representation of the true state PDF for a given problem It is in this context that the filter
PAGE 17
6 is used for the current research. While many of the above filters were shown to have perf orrnance superior to that of the EKF, realistic implementation may be difficult. Other techniques that have been proposed to solve the bearingsonly guidance problem include the modified gain pseudomeasurement filter (MGPMF)8 and modified gain EKF (MGEKF).9 Both use pseudomeasurements to avoid the need for linearization around a ''nominal'' trajectory, such as is found in the EKF. These methods are restricted to a special class of nonlinear functions, which can be manipulated into a pseudomeasurement form. The construction of implementable nonlinear filters for nonlinear stochastic systems remains a challenge. Numerous filtering algorithms, based upon series expansions to approximately realize the conditional mean, have been suggested (see the bibliography of Bucy et al. 10 ) Once again, these techniques represented a heavy computational burden, even for loworder dynamical systems. The SDREF may alleviate some of the problems discussed so far. The technique is derived from a new nonlinear regt1lator design methodology known as the state dependent Ricatti equation (SDRE) control method.1112 The objective is to take the nonlinear system described by x = f (x) + g(x)u and transform it to the statedependent coefficient form x = A(x)x + B(x)u where f (x) = A(x)x B(x) = g(x) (2 1) (22) (23)
PAGE 18
7 using direct parameterization. Once in this f ortn the SD REF is obtained by considering the statedependent coefficients to be fixed and applying the linear contint1ous Kalman filter to the ''frozen'' problem. The result is an algebraic SDRE whose structure is the same as the steadystate Ricatti equation. The emphasis of this research is to develop and analyze a guidance filter based on this technique. Initial work in the development of the SDREF1 3 demonstrates the filter in a simple two dimensional pendulum problem. Other application s include induction machine estimation, 14 simultaneous s tate and paramet e r estimation, 1 5 and gyroless satellite angula1 rate estimation. 16 The current research extends their work to include an exhaustive investigation into the internal properties of this filter in a strategic interceptor guidance problem as well as in the friction pendulum problem. Insight into the internal representation of the pdf when using parameterization versus linearization is examined, as well as the necessary conditions for stability. The first example presented here uses only nonlinear measurements, while the pendulum problem discussed herein adds the complexity of nonlinear dynamics as well. A third example demonstrate s the unique parameteriz a tion property of the SDREF when confronted with a loss of algorithmic observability
PAGE 19
CHAPTER3 SCOPE OF WORK This chapter outlines the scope of the work in this di sse rtation. Throughout the design and analysis of the SD REF, the EKF was used to compare perfonnance issues between linearized and parameterized methods Similarly a sampleddata bootstrap filter was used to compare probability distribution functions between the different nonlinear filters in the exoatmospheric guidance problem The three filter s are theoretically derived and applicationspecific equations are formulated for each problem. Controllability an d observability requirements for the SDREF algorithm are analyzed and determined Several performance measures are inve s tigated that include sensitivity to initial errors, noise parameters, and target maneuve1s. In addition, a somewhat nonstandard methodology, the bootstrap filter, is used to investigate the generally unobservable evolution of the pdf within the fi l ter The pdf was examined at distinct time steps to deter1nine whether the parameterized filter pdf is different from the linearized filter. All three filter s were coded in FORTRAN and installed in a 6DOF simulation for evaluation in an exoatmospheric terminal homing engagement. To examine the effects of both dynamic and measurement nonlinearities the SDREF and EKF were also coded in VisSim (a visual simulation tool for a simple, but highly nonlinear, pendulum problem ) Perfor1nance was measured in a manner s imilar to the exoatmospheric guidance problem. Because no f or1nal proof for stability could be found for the SDREF, modifications were made that allowed a detailed proof of local 8
PAGE 20
9 asymptotic stability to be found. The MSDREF is fully derived and an examination of its performance against the EKF for a problem consisting of loss of algorithmic observability is also explored.
PAGE 21
CHAPTER4 EXOATMOSPHERICGUIDANCEPROBLEM ENGAGEMENT SCENARIO The filters were initially analyzed in the context of an exoatmospheric ter1ninal homing engagement of an interceptor tracking a thrusting ICBM for approximately the la st 10 seconds before impact The interceptor was closing on the targe t at approximately 10 kilometers per seco nd. The target may perform thrust cutoff or other evasive maneuvers, which provide for a stressi ng nonlinear engagement scenario in which to test the filters. At the start of the engagement, the guidance filter was given erroneous initial state and covariance estimates. The sensor is assumed to have acquired the target and the guidance fi lter i s ready to estimate the relative motion between the two vehicles using passiveonly measurements. The estimates were u sed by a standard augmented proportionalnavigation ( PRO NAV) guidance law to command accelerations that keep the inte1ceptor on a collision course with the target. No errors associated with the attitude control sys tem were used as this would complicate the analysis and cloud the performance comparison. Figure 1 shows a model of the interceptor, a generic 4kilogram vehicle configured with divert and attitude control thrusters. A generic, SS18like ICBM model was used to generate the target motion. The interceptor may have many error sources, some of which include hand over errors, pointing errors, thrust errors, centerofgravity location errors, and seeker noise. Initially many of these noise sources were not used so that the filter characteristics could be more easily discerned. Several engagement closing 10
PAGE 22
11 angles were examined to allow for more stressing conditions (i.e., tail chase, noseon, beam shot). Once proper filter operation was assured, noise sources were added to more realistically emulate the true intercept scenario, as well as to identify weaknesses in filter performance. y Divert Thruster s Top View z Divert Thrusters Ba ck View y Attitude Control Thru s ter s X z z Figure 1. Generic intercept vehicle. The simulation used was the Guided Simulation Program Exoatmospheric (GSP EX0).17 It is a FORTRAN simulation developed by the General Electric Company for the Air Force. Figure 2 shows a block diagram representation of the subsystems modeled.
PAGE 23
12 Estimator _., Guidance LawN1 Azimuth Elevation Range Navigation State ~Commanded Estim:ue Acceleration Azimuth EJevatlon Range OwnState l o formation IMU AutoPilot Mixed Thruster Commands Vehicle Dynainics Figure 2. Guidance syste m block diagram To allow for rapid and efficient evaluation of filtering concepts, the author has modified the guidance system. It ha s a Monte Carlo capability with an outer loop around the Monte Carlo loop, to allow for parameter trades. Analysis was perfor1ned u s ing a personal computer~ however, analysis can be accomplished on other platforms, including the SUN, VAX, and CRAY.
PAGE 24
CHAPTERS EXOATMOSPHERICGUIDANCEPROBLEM FILTER DEVELOPMENT The following paragraph s outline the formulations of the EKF, bootstrap filter SDREF, and MSDREF that were used for thi s research. A brief theoretical des c ription of each filter is given, followed by the applicationspecific equations. There are three filter types tinder consideration in this research paper. The mathematical derivations for all of the estimators are fully developed and all assumptions regarding noise characteristics and linearization are explained. The filter s are presented in coordination with the problems that were used in examining their perfor1nance. The SDREF theory and derivation u se d in the exoatmospheric problem is shown in thi s Chapter the SDREF pendulum problem and derivation is presented in Chapter 7, while the MSDREF s tability proof and analysis i s found in Chapter 8. Many variations of SDREF are poss ible due to the infinite number of parameterization s available. Initially simple parameterization s were used. Another important factor in evaluating performance differences between the three filters was the selection of the noise characteristics. The EKF perf or111ance may be enhanced in some cases by arti s tically ''tuning'' the noise parameters within the filter To assure an unbia sed comparison between the filter s, noise parameter s were set to the same nominal values for all filter s. Due to several limitations on the bootstrap filter developed in a the s i s by this author, it was used only for pdf comparisons in the analy s is of the strategic guidance problem. The bootstrap filter needs further development to be useful in stressing 13
PAGE 25
14 engagements as the number of samples required for adequate pe1formance makes it intractable to run for large numbers of comparison. The SDREF and EKF were analyzed in great detail to assess key pe1formance characteristics including convergence, stability, and state estimation error. A key performance critique was the comparison of the evolution of the SD REF and EKF pdf' s to that obtained by using the bootstrap filter. Interesting insight into the differences between how the SDREF and EKF model nonlinearities was observed. There are many issues involved with the SDREF development, including correct parameterization of the nonlinear measurement equations, investigation into the requirements for solution of the SD REF, and the stability of the SD REF. These areas were explored to determine an accurate set of state and measurement equations and a feasible method to solve them. An initial non stressing scenario was selected to establish a benign comparative baseline performance assessment. The pdf' s were compared using tl1is baseline scenario as well as a nonlinear maneuvering target scenario. Once the baseline was established, the filters were tested against more stressing engagements and intercept angles. Performance was measured by examining the state estimate error and 3sigma standard deviations, as well as comparing the state estimate to the true state. The author has found that this is an accurate method of evaluating filter perfotmance. 5.1 EKF Derivation The derivation for the EKF can be found in many sources, 1819 and is presented here for completeness. The dynamics for the intercept problem at hand are linear, while
PAGE 26
15 the measurement equations are nonlinear. For this reason a combination of linear and EKF equations were used and implemented in discrete form. Given a linear timevaryin g differential equation for a dynamic system . x(t) = F(t)x(t ) + G(t)u(t) + L(t)w(t ) (5. 1 ) where x(O ) i s given, u ( t ) i s a determini stic input, and w(t) is a random di s turb a n ce, the corresponding discrete time s olution i s t k x(tk ) = (tt tki )x(tk 1 ) + J ( ti:, r) [ G (r)u( r) + L ( i)w( r)]di. ( 5 .2) t k1 If u and w are assumed piecewise constant over the sampling interval, this can be written as ( 5 .3) I t where r k t u k 1 = f (tk 'r) G (r)u dr ( 5 .4) t 1 1 and ( 5 .5) However if F, G and Lare also considered constant over the sampling interval, !lT, then Assuming then th at r(!lt) = (!lt)[I,, 1 (!lt)JF1G A (!lt) = (!lt)[/11 1 (!lt)]F1 L ( 5 .6) ( 5 .7) (5.8)
PAGE 27
16 the state propagation equation becomes ( 5 9 ) The covariance matrix can be propagated in a similar manner by ( 5.10 ) with ( 5 .11) To relate the covariance matrix of random sequences to the spectral den s ity of the random proce ss for the continuous disturbance input w ( t ) we note that E { w(t) w r (r) } = Q'c (t) 8 (ti) ( 5 12 ) The equation for Qkt then becomes tk Q k t = f ( tk, i) L ( i) Q 'c (i) L r (r)7 ( t k i)di( 5.13 ) 'k1 This now repre sents the integrated effect of the continuous disturbance input w ( t ) Thu s the two equations nece ss ary for the propagation of the s tate and covarian c e are ( 5 14 ) ( 5.15 ) where k t is defined by eFt:., and rkt is define in equation ( 5.7) The measurement process for the filter is now determined. Given the continuous nonlinear meas urement equation z( t ) = h(x(t ) ) + v(t) ( 5 16 ) The discrete version can be written as ( 5.17 )
PAGE 28
17 A complete derivation of the EKF update equations can be found in 18 however the final solution for the calculation of the Kalman gain, and the update of the covariance and s tate are where xk () and Pk ()come propagation equations ( 5.14 ) and (5.15) and (5.21) with E{vk} = 0 and E{vkvkr} = Rk. Now ( 5.14 ) and ( 5.15) and (5.18) through ( 5.20 ) define the full set of EKF equations. For the exoatmospheric intercept problem th e sta te consists of SXR syR szR VXR xk = VYR ( 5.22 ) VZR axT ayT azr where the states are relative position, relative velocity, and target acceleration The system equations are (5.23) ( 5 .24) (5.25)
PAGE 29
18 Putting thi s in state space notation, 0 SR 0 VR I 0 0 0 0 0 VR + 0 I 0 0 0 0 0 0 acom + 0 0 0 0 0 aT 0 0 0 0 0 0 I W T From this equation the s tate tran s ition matrix i s computed as I /lit (/it) = 0 I 0 0 where I is the 3x3 identity matrix I (eAt:.1 + Ailt 1) / .,1,2 / ( 1 e Mi) I ;t l eAtll (5.2 6 ) ( 5 .2 7 ) A ssuming the commanded acceleration remain s constant over the filtering int e rval J (tk, r)G( i)udi= The state propagation equation become s 1 2 flt AC0/11 I 2 /itAcon1J 0 ( 5.28 ) ( 5 29) The covariance propagation proceeds in a similar fashion The covariance can then be propagated by ( 5.30) The measurement equation is ( 5.31 )
PAGE 30
where and Range h k (x(tk )) = zE I Range YR I Range 19 i s the Gauss ian ( normal) random vector with zero mean and covariance R k (5.3 2 ) ( 5.33 ) This completes the derivation and equation form ul ation for the nines tate EKF needed for the interceptor guidance system. Given an initial estimate of the s tate x(O) and the covariance matrix P ( O ) equations ( 5.14 ) and ( 5 .15) and ( 5.18 ) through ( 5 .2 0 ) can be solved in a recursive fashion to achiev e an estimate of the relative vehicle motion. 5.2 Bootstrap Estimator Derivation The boots trap filter is derived u sing the recursive Bayes ian estimation relation s hip s. To begin the di scre te dynami c s tat e e quation is g iven by ( 5 .34) where fki is the system transition function and w k I i s a random, zero mean, white noise sequence representing a known calcu l ab l e probability distribution func tion At di sc r e te time s, an observation becomes available ( 5 .3 5 ) where v k is a zero mean, whitenoise sequence representing a known distribution function It is further ass umed that an initial pdf p(x0lfo) of the s tate vector i s available. The objective of this process i s to construct the pdf of the cu1Tent s tate xk given all the available inforrnation ; i e., p(xkl )
PAGE 31
20 The prior distribution before the update is defined as ( 5.36 ) where the probabilisti c model of th e sta t e evolution can be obtained throu g h the use of the syste m s tate e quation s and the known s tati s tic s o f wk as ( 5 .37) Because the proc ess noise i s a ssumed to be independent of any previous s t a t e valu es, deterministically The update stage occur s when a measurement becomes available The a priori den s ity may be updated according to B ayes' rule as (5.38) where the normalizing denominator i s (5.39) and the conditional pdf, p (yk lxk), is defin e d by the measurement model and the known statistics of v ,,_ The method by which the fallowing th eory i s implemented in the bootstrap filter will now be defined The boots trap filt e rin g algorithm provides the mechani s m for defining a set of random s ample s {xk_1 (i); i = 1 ... N} from the pdf and obtaining th e di s tribution p(xk_1 I_ ,) .5 7 The propagation stage assumes a number of N random sample s have been tak e n from so m e initial pdf p (xk _1 I,). Each random s ample for
PAGE 32
21 each state is then passed through the system model to obtain a set of a priori s ample s at t ime step k x k (i) = fk ( x k J (i), w k I ( i )), i = 1 ... N ( 5.40 ) where the a s teri s k repre s ent s the a priori sample. The update s tage of the filter take s each of the a priori samples and compute s a se t of a po s teriori s amples from which the estimate can be obtained. To do thi s, the likelihood of each a priori s ample l ( x k I yk) i s calculated using the measurement equation s and a normalized weight is given to each sample by ( 5 .41) This defines a dis c rete distribution over the a priori den s ity. The normalized value s of q ; can be resampled u s ing a uniform di s tribution. A larger value of q ; allows more sample s from that bin to fall into the a posteriori pdf. The state associated with each new sample selected is assigned to the new a posteriori state {xk ( i); i = 1, . N}. Thi s a po s teriori s et of samples is approximately distributed as p ( x k I) From this the mean and variance for each state estimate can be computed. Because this filter makes no assumption of linearity or Gaussian statistics, it was useful in providing a close approximation of the true pdf Thi s s emitrue pdf can then be u s ed to inve s tigate the effect of linearization a nd parameterization on the evolution of the pd f
PAGE 33
22 5.3 SDREF Derivation The SDREF is the nonlin ear ''dual' of a nonlinear regulator control design tec hniqu e which u s es a parameterization technique to bring the nonlinear dynami c and measurement equations X = j (x) z = h (x) and tran sfor1n them into a linear like s tatedependent coefficient (SDC ) form x= F(x)x z = H(x)x ( 5 42 ) ( 5 43 ) ( 5 .44) ( 5 45 ) The distinctive ass umpt ion in the SDREF technique is that the SDC s can be treated a s ' frozen'' ove r a smal l interval to give adequate piecewise performance. The algebra ic Ricatti equat io n is s o l ved at every fi l ter measurement update and t h e linear Kalm an filt ering algori thm s applied. The technique was first applied to the fol l owing nonlinear regulator problem .13 The ta s k is to minimize 1 00 I= 2 J x7 Q (x)x + u7 R(x)u dt to ( 5.46 ) wit h re s pect to the s tate x and contro l u, subject to the nonlinear differential constraint i = f (x) + g(x)u Th e SDREF approach for obtaining a sub optimal solution i s t o use direct parameterization to bri n g the nonlinear dynamic equat ion to the SDC form x = A(x)x + B(x)u ( 5 47 ) ( 5 48 )
PAGE 34
23 where f (x) = A(x)x and B(x) = g(x) From this point the SDREF can be solved for P as A r (x)P + PA(x)PB(x)R 1 (x)B r (x)P + Q(x) = 0 ( 5 49 ) and the nonlinear feedback controller can be co nstructed as u = R1 (x)B r (x)P(x)x (5.50) The SD REF technique is shown 1 3 to have desirable properties of local asymptotic s tability and suboptimality given mild conditions of stabilizability and detectability. It also shows that if A1 (x) and Ai (x) are distinct SDC parameterizations they can be combined to yield a third parameterization by A(x, a)= a41 (x) + (1a)Ai (x) This technique of co mbining parameterizations allows an infinite number of parameterization possibilities to enhance controller design. (5.5 1 ) An example1 3 of the SDC form and SDREF is repeated to introduce the method. Consider the stochastic nonlinear system i=f(x)+f(w) z = h(x) +v ( 5 52 ) (5. 53 ) where w is Gat1ssian zeromean white process noise with E[w(t)wr (t + r) ] = Q (t)8(1:) and v i s Gaus s ian zeromean white measurement noise withE[ v(t)vr (t + r) ] = R(t)o( i). After bringing the system to the SDC form i = F(x)x+ fw z = H(x)x+ v ( 5 54) (5.5 5 )
PAGE 35
24 t h e SDREF i s given by the s teady state, linear continuou s Kalman filter equation s a s x = F(x)x+ K_(x)[z(x)H(x)x] where and P i s the positive definite solution to the SDREF F(x)P + ppT (x)PHT (x)R 1 H(x)P + r r Qr= 0 (5.56) (5.57) (5.58) The preceding method was applied to a s imple twodimen s ional pendulum problem .13 Res ult s were comparab l e to those from an EKF. The main thru s t of the re searc h was to develop a SDREF for the linear dynamic and nonlinear measur e ment equation s associated with the strateg i c, terminal homin g gu idance problem and i nve s tigate filter charac teristic s u s ing the EKF and bootstrap fi l ter for comparison and ana l ysis. 5.3.1 AngleOnly Measurement SD REF The following derivations are perform e d for the bearing measurement s only SDREF. It shows that thi s form of the SDREF i s not d e tectable and ther efo r e viol a te s the neces s ary conditions to assure the positive definite solution of the Ricatti equation. Both the continuous and discrete forms of the filter are given, as well as s tabilizability and detectability analy ses for each type The filter derivation begin s with the definition of the s tate e quation s .x = Fx +Gu+ L w (Co ntinuou s form ) (5.59) (5.60)
PAGE 36
and the measurement equations The states to be e st imated are 25 z = h(x) + v (Non linear form) z = H x + v ( Linear form ) SXR syR SZR VXR xk = VYR VZR axr ayT a zr where the s tate s ar e relative position, relative velocity, and target acceleration. The system equation s are V . Ar A R ron1 Putting this in state space n otation 0 I 0 SR SR 0 0 0 0 I VR VR + I aco, n + 0 W T 0 0 }J aT aT 0 I X [ F] {x}+ [G] u + [L] w The measurement equation i s (5.61) ( 5 .62) ( 5 63 ) (5.64) (5. 65 ) ( 5.66 ) ( 5 .67) ( 5 68)
PAGE 37
26 where, based on the smallang l e approximation, h(x) = and zR Range YR Range The SDREF approach for obtaining a sub optimal s olution i s to use d i rect parameterization to bring the nonlinear dynamic equation to the SDC form where and i = F(x)x + G(x)u f (x) = F(x)x z = H(x)x+v From thi s point the SDREF can be solved (5.69) (5.7 0 ) (5.7 1 ) F(x)P + PF7 (x)PH7 (x)R 1 (x)H(x)P + LQL7 (x) = 0 (co ntinuou s) (5.72) Since the system dynamic equations are linear the only parameterization needed i s for the measurement equation. From the definition z= Hx+v the linear form for the SDREF can be cast as H = 0 0 1/ R O O O O O 0 0 1 / R O O O O O O 0 5.3.1.1 C ontinuo u s SDRE F Using the abovedefined s tate and measurement equations with A= 0 1 the following matrices are used for the filter equations ( 5 .74) (5.7 5 )
PAGE 38
27 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 F= 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 .1 0 0 0 0 0 0 0 0 0 .1 0 0 0 0 0 0 0 0 0 .1 The time constant for the Markov process A, was base d on an expected target acceleration of approximately 981 (m/ s2 ) 0 0 0 0 0 0 0 0 0 1 0 0 G = 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 (5.76) (5.77) The proce ss noise u se d for this problem is 1 g on the xtarget acceleration and 3 g's on both y and z target accelerations. The disturbance input is E[ w(t)] = 0 E[ w(t), w7 (r) = Qc(t)8(tr) ] ( 5 .78) ( 5.79)
PAGE 39
28 which for small values of 11t can be expressed as 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Q k 1 = LQ'c L r 11t = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 11t ( 5.80 ) 0 0 0 96 0 0 0 866 0 0 0 866 where 866 ( m I s2 ) is the variance of 3 g's of acceleration noise in they and z direction s The parameterized measurement matrix can be computed as H(x) = 0 0 0 1 Ran ge 1 0 0 0 0 0 0 Rang e 0 0 0 0 0 0 0 The standard linear continuou s Kalman filtering equations follow as = F(x)x + K1 () [ z(x)H(x)x] where and P i s the pos itive definite so lution to the s tatedependent Ricatti equation (5.8 1 ) (5.82) ( 5 .83) where Q k i = LQc L7 !:lt. Ricatti equation theory provides for a legitimate pos itive semi definit e so lution if the following co ndition s are m e t 0 R = R T > 0 m < n, X E n ( 5. 8 5 ) and (F, H) is pointwise detectable in the linear sense and( C F ) i s point wise
PAGE 40
29 stabilizable in the linear sense for all x E Q where C is a fullrank factorization (FRF ) of Q k 1 given by ( 5 86 ) For the given linear dynamics we can construct the controllability matrix as Con= [ C FC ... p n 1c] ( 5 87 ) If Con is of full rank, controllability and therefore stabilizability is assured. If it is not fully controllable, a minimum condition of stabilizability must be met. Putting the s ystem equations into the controllability canonical form the uncontrollable poles c a n be checked for s tability 20 Using a similarity transformation, the controllability canonical fonn can be found as F = TFT C = TC where Tis unitary and the transformed system has a staircase form F,,c F = 0 0 (5 88 ) (5 89 ) where ( F e G e ) are controllable and ( F11c ) i s uncontrollable. The poles of the uncontrollable subspaces are then checked to insure stability. If they are found to be stable, then stabilizability i s verified. Using the noi s e matrix in equation (5.13) and performing the FRF for C it can be shown that the controllability matrix bas full rank and the ref ore stabilizability is verified.
PAGE 41
30 Similarly, the observability Gramrnian matrix is constructed as H HF Obs= (5.90) HF" 1 If Obs is of full rank for all x, then full pointwise observability in then linear sense is assumed. Once again pointwise detectability needs to be assumed for the Ricatti equation solution to be positive semi definite. If it is not of full rank the similar observability canonical form can be used to again check the pointwise stability of the unobservable poles. This canonical form is given by F;,o F= 0 F 2 1 [ H=O F 0 (5.91) Using the above referenced matrices for F and H, it can be shown that the observability matrix is not of full rank. Therefore, pointwise stability of the unobservable poles must be checked. It was found that the poles were not pointwise stable and therefore point wise detectability was not verified. A similar analysis was then perfor1ned on the discrete form of the filter. From this point on, only the discrete form of the SDREF is analyzed. This for1n was chosen due to the use of discrete measurements that rue only available at a subset of the integration step size of the simulation. The discrete filter is better suited for this type of simulation and would most likely be used in any realistic guidance software.
PAGE 42
31 5.3.1. 2 D i sc r e t e SDREF The discrete form of th e SDREF i s d e riv e d u s in g the follow ing formulas. Since F i s a constant matrix the s tate tran s ition matri x can be found t1sing Laplace tran s form s as I Jilt (llt) = 0 I 0 0 l(eJ.i>J +Mtl)I 22 J ( 1 e Mr) I /4 leAlli When u s ing 2 = 0 1 and flt= 0 025 the following tran s ition matrix is computed 1 0 0 .025 0 0 000312 0 0 0 1 0 0 0 2 5 0 0 000312 0 0 0 I 0 0 0 2 5 0 0 .000312 0 0 0 1 0 0 .0249 0 0 = 0 0 0 0 1 0 0 0249 0 0 0 0 0 0 I 0 0 0249 0 0 0 0 0 0 .997 5 0 0 0 0 0 0 0 0 0 .99 75 0 0 0 0 0 0 0 0 0 .997 5 As s uming the commanded acceleration remain s constant over the filtering interval ruk i = J (tk, r)G(r)ud"!' = 1 2 f:J Aeon,/ 2 fltAcun,/ The s tate propagation equation becomes 1 2 flt I 2 0 xk = [ ]xk1 + flt! Acom 0 (5.92) ( 5.93) ( 5 .94) (5.9 5)
PAGE 43
32 Also, (5.96) and t (5.97) Using 1 g of acceleration noise on the x target acceleration and 3 g's in y and z and performing the above integration the process noise matrix is obtained. 9.46e9 0 0 9.36e7 0 0 4.99e5 0 0 0 8.54e8 0 0 8.44e6 0 0 4.45e4 0 0 0 8.54e8 0 0 8.44e6 0 0 4.45e4 9.36e7 0 0 9.98e5 0 0 .006 0 0 Qk1 = 0 8.44e6 0 0 9.0e4 0 0 .054 0 0 0 8.44e6 0 0 9.0e4 0 0 054 4.98e5 0 0 .006 0 0 .478 0 0 0 4.Se4 0 0 .054 0 0 4.32 0 0 0 4.5e4 0 0 .054 0 0 4.32 The discrete Ricatti equation solver requires the following conditions to be met to in s ure that a legitimate positive semi definite solution to the Ricatti equation is obtained: T' Q, p E R'IXl1 H T E R',xm' RE Rmxni' Q =QT> 0, R =RT> O,m < n, XE n (5.98) and (,H) is pointwise detectable in the linear sense, and (C,) is pointwise stabilizable in the linear sense for all x E Q, where C is a FRF of Qki and is defined as ( er C) = Qkl and rank( C) = rank(Qk t) (5.99) It is also must assumed that is invertible without need for further modifications to the solver. If all assumptions are met, then the Ricatti equation will have a unique non negative definite solution throughout Q
PAGE 44
33 Using the above referenced matrices, the same analysis used for pointwise stabilizability of the continuous system was perforrr1ed. The controllability matrix of the pair ( C, ) was found to be of full rank and therefore the system is stabilizable. Again an analysis similar to the continuous version was perforrned using the pair (, H) and the pointwise observability matrix was again found to be not of full rank The canonical form was calculated and again the poles of the unobservable st1bspace were determined to be pointwise unstable and therefore the system cou ld not be declared pointwise detectable in the linear sense. Because the system has been shown to be pointwise undetectable, the Ricatti equation solver cannot provide the positive semi definite solution to the Ricatti equation. Consequently, another app ro ach to the problem must be examined. One way to make the system detectable is to allow a measurement in the range or xdirection This is shown in the next section. 5.3.2 SDREF With Range Measurement This section analyzes the same filter described above except that a range measurement has been added to the measurement equation. By adding the range measurement the system becomes both stabilizable and pointwise detectable under the conditions set forth above. Two types of range measurements were evaluated to provide adequate observability for the system.
PAGE 45
34 5.3.2.1 Intensity range measurement The range measurement used is the so urce radiant inten si ty of the target arriving at the focal plane and is where lk = intensity (watts) J;,, = source radiant intensity ~A = bandwidth (m) watts sr m r = radius of seeker aperature (m) R = range from source to seeker (m) v k = white Gaussian noise This gives the measurement equation: h(x) = z / Range yl R ange J;., ~:im2 I Range2 Putting equation (5 101) into the SDC form 0 0 j 2 2 2 Range =xR + YR + ZR 1 0 0 0 0 0 Range 0 H= 0 1 Range 0 0 0 0 0 0 0 l,i~/4m2 0 0 0 0 0 0 0 0 Range2 xR (5 100 ) (5.101 ) (5.102) The measurement noise for the range measurement is somewhat more complicated as it is range dependent. The actual intensity measurement has Gaussian noise from the focal plane electronics, which can cause error in the measurement. However, it is assumed this error is small in comparison to the fluctuations of the source radiant intensity and is
PAGE 46
35 assumed to be zero. l;., is however comprised of both the deterministic source radiant intensity and a random part given by (5.1 0 3) The resulting inten s ity contribution from the random fluctuation s in the plum e core can be calculated from equation ( 5. 101 ) as I = J;., ran ~},. :n:r 2 ran Ran ge2 ( 5 104 ) T he variance of this random component of intensity is obtained by taking the expected value of the square of Iran (i.e., E{Jran 2 } ) The equation for the expectation operation gives E 2 l,i LlAm2 ran Ran ge2 ( 5.105 ) If p = E { 12 ..i,m, } then the range depend ent measurement variance term for the filter becomes R K 2 p~m2 Ran ge2 ( 5 106 ) Using appropriate values for these constants, a meas urement matrix can be given by 2.5e9 R k = 0 0 0 0 2.5e9 0 ( 5.10 7) 0 8.9e7 Using the abovedetermined matrices, the SDREF was again calculated. The filter responded well during the first three to four seconds of flight and then diverged The
PAGE 47
36 problem lies in the instabi l ity of the Ricatti equation solution. In examining the point wise c l osedloop eigenval u es of the system (i e., the spectrum of (5.108) with x, in this case range frozen), the system would be asymptotically pointwi se stable in the linear sense if the closedloop eigenvalues all lie inside the unit circle. The eigenvalues of the SDREF were computed for each filter update stage. They began near the unit circle (i.e., 0.9992) and grew to become unstable as the s imulation ran. The final step in the Ricatti equation solution is the solution of an nth order linear matrix equation. A number was computed to determine the condition number for the inverse matrix. It could be seen that the condition number increased (i.e., illconditioning) as the closed loop eigenvalues approached the value of one. The relatively small size of the range variance causes R 1 to become very large, and therefore il l conditioned. Once the closed loop eigenvalues became larger than one, the problem diverged quickly. Because of this numerical instability a different range measurement was tried. 5. 3.2 2 D is tan ce ra n ge m eas ur e m e n t The new range measurement is the distance measurement in meters. This value is assumed to be available from the signal processing algorithms. The new measurement is (5.109)
PAGE 48
37 Therefore, the new measurement equation becomes z Range h(x) = y Range Rang e and the measurement matrix becomes 0 0 1 Range H= 0 1 0 Range Range 0 0 XR (5.110) 0 0 0 0 0 0 0 0 0 0 0 0 (5.111) 0 0 0 0 0 0 Using these values in the SDREF, satisfactory results were obtained. (These re sul t s are given in Chapter 7.) It was decided that thi s would be the formulation that would be used to further test this new filtering technique. Before the se results are presented, however, an investigation into the theoretical st ability of the SDREF is presented. 5 .3.3 T h eo r e ti c al Inves ti g a ti on s One of the main theoretical criteria for filtering techniques is the stability of the algorithm. The fallowing section reviews several theoretical stability inve s tigations. While these investigations were unsuccessful in proving stability, it is hoped that they provide valuable insight into the SDREF. It i s these attempts that led to the development of the MSDREF s hown Chapter 9. While the methods outlined below approach the problem from seve ral different theoretical aspects, all were unsuccessful in s howin g local stabi l ity for the SDREF. Because the true and estimated s tates x and x both appear in the error eqt1ation
PAGE 49
38 e = (F(x)xF(x)xK(x)[ H(x)xH(x)x]) (5,112) each of the approaches fail when presented with the cross tertns that are created. It i s hoped that this section will prove useful to future stability investigations. The approaches are shown s tarting with s imple linearization techniques and progress through more complicated and sometimes restrictive theorem s. The first of these is a s imple Taylor series expansion of the parameterized equations 5.3.3.1 Taylor series expansion The following stability analysis will exploit a Taylor series expansion about the local zero point. Consider the stochastic, nonlinear system x = f (x) + g(w) z=h(x)+v (5.113) ( 5.114 ) where w i s Gaussian zeromean white process noise with E[ w(t)wr ( t + i)] = Q(t)8 ( i) and v is Gaussian zeromean white measurement noise with E[ v(t)vr (t + i) ] = R(t)8 ( r). After bringing the sys tem to the SDC form x= F(x)x+rw z = H(x)x+v ( 5.115 ) ( 5.116 ) the SDREF is given by the linear continuous Kalman filter equations having s tate dependent coefficients x = F(x)x + K(x)[z(x)H(x)x] where K(x) = P(x)H r (x)R 1 and Pis the positive semi definite so lution to the algebraic SDRE F(x)P(x) + P(x)F r (x)P(x)H r (x)R 1H(x)P(x) + r r Qr= o ( 5.117) (5. 1 18) (5.1 19 )
PAGE 50
39 Conjecture: Assume that R > 0, Q > 0, f (x), and h(x) ERK, and that the parameterizations F(x) and H (x)can be shown to satisfy that ( F(x), H(x) )is a pointwise detectable pair in the linear sense and ( C(x), F (x)) is a pointwise stabilizab l e pair in the linear sense, where CC T = Q k i. Also assume F(x) E C00 and H(x) E C00 are analytic where C00 is the class of functions whose derivatives of all orders are continuou s, and that f (0) = h(O) = 0. Then the SDREF is locally asymptotically stab le. Outline of Attempted Proof: For ease of notation the proof will be examined in scalar form. Given that the scalar system dynamic equations are x= F(x)x+Gu+Lw z = H(x)x+v and i = F(x)x +Gu+ K[ zH(x)x] ( 5 .12 0 ) and defining the estimat ion error rate as = .x x the above equations can be written as E = (F(x)xF(x)xK(x)[H(x)xH(x)x]) ( 5 .121) Now expanding h(x) and H(x)x in a Taylor series expansion about zero h(x) = H(x)x (5.122) ciJ1(0) a2 h (xo)2 h(0)+1 (x0)+ 2 +H. 0 .T.= dx x=O ax, x=O 2 8H(O) a2 H H(O) +1 (x0) + 2 ax x=O ax x=O 2 ( x o) + H 0 T x 2! Also for f (x) = F(x)x f (x) = F(x)x ( 5.123 ) q'(O) J2J ( x o )2 f(0)+1 (x0)+ 2 +H.O.T.= ax x=O ax x=O 2 8F(O) J2 F (xo)2 F(0)+1 ( xo ) + 2 +H.0. T x ax x=O dx x=O 2 S i milar equations can be written for h(x) = H(x)x and f (x) = F(x)x.
PAGE 51
40 Equating like terms Jh(O) = H (O) = Jh(O) dX x=O dX x=O and J/(0) = F(O) = J/(0) dX x=O dX i=O Using this and recalling that K = PHT R 1 e(t) = (F(O)PHr (O)R 1H(O)))e(t) + H.O.T. (5.124) Herein lies the problem, the higher order terms contain both x and x. These terms cannot be guaranteed to approach zero, much le ss go to zero faster than the first ter1n. If the error dot equation could be expressed in this linear form a local stability proof could be achieved in the following manner. For ease of notation, let F(O) = F, H(O) = H, and B(t) = Following a similar derivation, 19 a Lyapunov function is chosen ( 5.125 ) where I is the positive definite information matrix, the inverse of P. From the definition of the derivative of the inverse, By substituting this into the matrix Riccati equation, a differential equation i =IF+ pT I+ ILWL T I HT N 1H is yielded that has the steadystate solution 0= IF+ F T I + ILWL T I HT N 1H ( 5.126 ) (5.127) (5. 128 )
PAGE 52
41 The time derivative of the Lyapunov function is V =2Tle=2T[FKH] = 2TI[FPHr N 1H]c = cT[(IFHr N 1H)+(IFHr N 1H)]c Using (5.128) this can be written as (5.129) (5.130 ) As the terrn in square brackets is positive definite at all times, V < 0 at all times. This would verify the local asymptotic stability of the SDREF for some small neighborhood about the origin. The estimation error would approach zero as time approaches infinity. 5.3.3.2 Using a Taylor series expansion about the true states This approach varies slightly from the above, in that a Taylor series expansion of the estimated dynamics and measurement equations are linearized about the true states. Using the same state and measurement equations (5.1135.120) and the error equation (5.121) in the preceding proof, F(x) and H(x) are expanded as F(x) = F(x)+F'(x)(xx)+ F''(x) (xx)2+ ... H.O.T. 2 A ) H( ) ,. ) H''(x) ) 2 T and H(x = x + H (x)(xx +(xx + ... H.O. 2 Using equations (5.131) and (5.132) equation (5.121) is written as = { F(x)KH(x)}c+{ F'(x)KH'(x) } 7J(2 ) (5.131) (5.132) (5.133) Once again the problem with this approach can be seen. The only way it can be shown that the second term disappears or is negative is for to be zero. This can only be shown to be true for the linear case.
PAGE 53
42 5.3.3.3 Taylor series expansion about the error Again, using the same system equations ( 5 113 5 121 ) as in the first Taylor series approach, this second approach linearizes about the error. Linearizing the dynamic and measurement equations F''(s) F (x) = F (s) + F'(s)(xs) +(x)2 + ... H 0 T 2 F'''(s) F'(x) = F'(s) + F''(s)(xs) +(xs)2+ . H.O. T. 2 H''(s) H (x) = H (s) + H'(s)(xs) +(x )2 + ... H 0. T 2 H'''(s) H'(x) = H'(s) + H''(s)(xs) +(xs)2+ . H O.T. 2 After s ub s tituting these into error equation ( 5 .12 1 ), e = (F KH)s + (F' KH')(xs)s + (F' KH')ci + (F'' KH'')(xs)ci + (F'' KH'')(x )2 2 + (F' '' KH''')(xs)2 x + ?J(s2 ) 2 (5 134) (5.135) (5. 1 36 ) (5.137) (5.138) Again the error equation i s only stable if it can be shown that 0. This can only be shown to be true for the linear case. 5.3.3.4 Residual error stability approach The followin g approach illu strates, under conditions of pointwi se detectability in the linear se n se th at ensure a legitimate so lution of the SDRE, th at the SDREF can be shown to b e asymptotically sta bl e .21 The system dynamic s and filter equations are given in discrete form as
PAGE 54
43 with ( 5.139 ) ( 5 140 ) ( 5 .141) ( 5 142 ) ( 5.143 ) ( 5 144 ) ( 5 145 ) where ~ +1 (xk+i) i s the solution of the SDRE. These are the equations necessary for the discrete SDREF. It is also assumed that vk+J and wk+i are random, zero mean inputs and therefore, the stability of the SD REF is independent of these error sources. The state estimation error both after and before the update is defined respectively as " X k + I = Xk + l X k+I " x k +llk = xk+1 xk+llk Also for ease of notation And finally a candidate Lyapunov function V k + i is defined as V. TP, 1k + l = xk+I k+I xk+ I ( 5 146 ) ( 5 147 ) ( 5 148 ) ( 5 149 ) ( 5 150 ) ( 5 .151) The conditions for which V k + i is a decreasing sequence must be determined therefor e showing the exponential stability of the SDREF.
PAGE 55
44 In the preceding stability approaches, linearization of the filter equations has not shown stability of the filter. This attempt constrains the state transition matrix and measurement noise matrix so as to characterize the error due to linearization. The measurement res idual error can be characterized as ek + I = Zk+I h(Xk+llk) ek+ I = H k + I (xk+I )xk+I Hk+I (xk+I )xk+I (5.1 52 ) (5. 153 ) Rewriting Hk+i (xk+i) using a Taylor series expansion, for each output error prediction component of ek+I namely ei,k+l' there always exists a residual 't,k+t due to the first order linearization in order to obtain an exact equality. Therefore e i,k+ I = H i,k+I (xk+I )xk+ I {Hi,k+I (xk+l) + H'i,k+I (xk+I )(xk+ I xk+1) + ... }xk+l ( 5.154 ) and the residual error is 't,k+I = {H;i,k+I (xk+I )(xk+I xk+ l )+ . } xk+I Then the exact measurement error can be defin e d as or written in another form where a i,k+I i s an unknown timevarying scalar related to r; ,k+t by 1i,k+1 = (1a i,k+l )ei,k+J Introducing the signal vector where ak+ I E R P P i s an unknown time varying diagonal matrix defined a s ak+I = diag(a;,k+ 1, . ,ap, k + 1) ( 5.155 ) (5.156) ( 5.157 ) ( 5 158) (5. 159 ) ( 5 160 )
PAGE 56
45 If both sides of equation ( 5.142 ) are subtracted from xk+l' P. H T R 1 xk+I = xk+Jlk k+t k+1 k+t ek+1 ( 5 161 ) Sub s tituting this into equation ( 5 151 ), and u s ing the fact that Pk+t = Pk+i T l l TP. 1T H T R 1 y k+I = Xk+llk k+l Xk+llk Xk+llk k+I k+1 ek+I T R IH T R 1 H P. H T R 1 ek+I k + l k+lxk+l/k + ek+I k+I k + l k+I k+I k+I ek+I ( 5 162 ) Because the c ovar i ance matrix i s not propagated during the time between updat es, it i s assumed that the a priori covariance before update i s equal to the previou s updated covariance ( 5 .163) Therefore P. I P, 1 k+I = k ( 5 16 4) Before proc e eding f urther th e inver s e relation s hip i s defined a s P. I P, 1 H T R I H k+1 = k+Ilk + k+I k+I k+I ( 5.165 ) Once again a s imilar problem ari ses Equat ion ( 5.165 ) i s not derivable f rom th e bas i c premise of the SDREF update equations .21 If thi s equation was available it could b e shown that Vk+i is a decrea s ing s equence. Then l l l l T H T R 1 v k+t = v k+tlk xk+tlk k+1 k+1 ek+1 T R l H T R I H P. H T R 1 ek+l k+I k+1Xk+llk + ek+l k+l k+l k+I k+1 k+l ek+ I ( 5.166 ) and re v iewing th a t ( 5 167 ) and T H T T T x k+1 k+t = e k+1a k+1 ( 5 168 )
PAGE 57
46 Equation (5.166) is written as ( 5 169 ) And using the fact that V T D TP, I p k +ltk = xk r k k+l/k kxk ( 5 170 ) with Equations (5.167168), Equation (5.169) is rewritten TT T T T { T R 1 R I R l r k+I r k = ek+t ak+t k + t ak+i ak+1 k + 1 k+1 ak+1 + Rk+i1 Hk+1 Pk+1Hk+1 r Rk+i1 }ek + t ( 5.170 ) T{r:;,TP,lF, P,1}+ xk r k k k k xk Therefore under the condition of pointwi se detectability in the linear sense, and certain restriction s on R and Q (i.e., R k > 0 and Q k 0 ) then, and the system is asymptotically stab l e a s the values in parentheses are 0. 5.3.3.5 Lyapunov approach using the expectation operator This stability approach uses a similar method.9 Given the system equations .x. l = L1 .x. + w t+ I I x = x + K [H(x )x. H(x. )x. ] + K v r I r I I I I II the error can be updated and propagated by the equations e ; = X ; X; (Updated error) e; = X; X; ( Propagated error) ( 5 .17 1 ) ( 5 .172) ( 5 .173) ( 5.174 )
PAGE 58
47 Using equations ( 5 171 5.172) e = x x K [H (x )x. H(x. )x. ]K v I I I I I I I I II e = A ,e I + w I I ,e.i.;_ I /( 5 175 ) (5. 176 ) After some manipulation of these equations the erro1 update equation can b e developed e.+1 = A e A K H1x + A K H2x K v + w I I I ,e~ I I I I I I I I where H1 = H(x) and If the following Lyapunov function is defined V Ti+I = e;+i e;+1 it can be shown that E{+1 (e;+I) (e;)} 0 Substituting equation ( 5.177 ) into equation ( 5 179 ), the following equation can be computed E{+1
PAGE 59
48 This same approach was also tried with the Lyupanov function V=er e and Equation (5.182) was developed E{V:+1 (ei+I )(e;)} = e;rE{ Ar Al}ei +E{2.xAT AKH1x2xr A r AKH1x} +E{2xT AT AKH 2x+2xT AT AKH 2 x} +E{2x7H,TKT AT AKH2 x} + E{xrH1r KT KH1x} +E{xTH 2 T KTKH2x} ( 5.181 ) (5.182) Once again cross terms are generated and cannot be shown to be negative definite. All of the previous attempts at stability proofs fail for the same reason: The cross ter1n xx, produced in the error equation, raises significant difficulties in finding a closed form differential equation for Attempts to limit the problem to a scalar were not able to resolve the cross coupling problem. While several standard and nonstandard approaches were tried, a satisfactory stability proof could not be verified It is hoped that these investigations can prove useful in future attempts to deter1nine filter s t ability.
PAGE 60
CHAPTER6 EXOATMOSPHERICGUIDANCEPROBLEM SIMULATION RESULTS Chapter 6 investigates the comparative performance of each filter against a variety of engagement conditions and noise parameters. Before any study on filter perfor1nance was done, a baseline perfor1nance analysis was completed, which allowed a detailed investigation of the SDREF against one engagement scenario. Once the baseline performance was established, further investigation into filter sensitivity, robustness, and closedloop performance was performed. A study of pdf evolution for each filter is also shown. 6.1 Base l ine Performance To establish a baseline performance measure, an engagement representing the final 10 seconds of a closing interceptor/boosting ICBM engagement was chosen. Table 1 details the major initial filter parameters. The aspect angle is defined as the initial Euler rotations between the axes of the interceptor and the target. For example, an aspect angle of (30,30,30) would signify that the interceptor's initial velocity vector has been rotated 30 degrees about each axis. An aspect angle of (0,0,0) would be a headon engagement. The (30,30,30) aspect angle was chosen so that all axes of the filter would be exercised in the engagement. SIGP and SIGVP are the standard deviations for the random initialization errors on the ''true'' position and velocity states respectively. When SIGP 49
PAGE 61
50 and SIGVP are zero, the filter has been initialized with the ''true'' initialintercept conditions. SIGXHA T are the values used to initialize the covariance matrix for the EKF and the bootstrap filter. XHAT is the initial state given to the filters. The measurement noise for the xdirection is a rangedependent noise that is given as a percentage of the range (i.e as the interceptor gets c loser to the target the noise on the range meas urement decrease s). All filters were made to operate with the same initial conditions and statistics so that comparisons could be made. A perfect filter was used to guide the vehicle to impact, while the EKF, SDREF, and bootstrap filter were run open loop. The bootstrap filter was run with 10,000 samples. Table 1. Filter simulation parameters baseline conditions Parameters ( units ) X y z Aspect Angle ( degrees ) 30 30 30 SIGP ( m) 0 0 0 SIGVP (m/s) 0 0 0 (J'2 (m2) l.OE5 1.0E5 1.0E5 r SIGXHAT (J' 2 V (m2 / s 2) 1.0E3 1.0E3 l.OE3 (J' A 2 (m2/s4) l.OE2 1.0E2 l .OE2 R (m) l.OE5 10.0 10 0 XHAT V (m/s) l .OE4 335 9 223.9 A (rn/ s2) 55.2 57.40 38.2 SIGQ (EKF) (m2 / s4) 196.0 196.0 196.0 SIGQ (SDREF) (m2 / s4) 1.0E6 l.OE3 1.0E3 Meas Noise (%,m,m) 0.1 50. 0 50.0
PAGE 62
51 The following graphs reflect filter performance. Figures 3 through 5 show the positionstate estimate errors with .fci standard deviations for each filter. Velocity and acceleration perfor1nance was similar between the EKF and SDREF. The bootstrap acceleration estimates are poor due to an inability to forward mea s urement inforrnation to the second derivative states. The problem can be remedied by the use of millions of samples or other adhoc improvement schemes; however computer limitations and the scope of this dissertation do not allow inve s tigation at this time Figures 3 through 5 give initial insight into the relative performance of each of the filters for the given scenario. All of the ftlters adequately estimate the state to the precision necessary for an intercept to occur With state estimate errors remaining within a standard deviations, the filter s perform well. The bootstrap filter does show a tendency to keep a wider variance than the other two filters. However as discussed in paragraph 6.2, when the EKF and SDREF are run in Monte Carlo fashion, they too have wider variances. The results shown here are for one run, while noticing that the bootstrap filter is inherently Monte Carlo even for one run. The performance curves in Figure s 3 through 5 verify proper filter operation and establish a baseline perfor111ance from which more indepth analysis proceeds It is intere s ting to note in Figures 6 and 7 that the standard deviations for the EKF and SDREF are approximately identical after the initial transients in the EKF die out. This indicates that once the EKF has settled and the estimates are being computed accurately, the EKF covariance propagation and update are equivalent to the algebrai c Ricatti solution being computed in the SDREF. This relationship is only shown for one state; however, similar results were found for all states.
PAGE 63
E g a. )( e g "r l N 200 100 0 100 1 0 0 10 0 200 100 0 )( 100 200 15 10 e 5 0 l .5 >10 15 15 10 s E 0 g ".5 N 1 0 200 100 .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. E 0 i 100 )( 15 1 0 5 E 0 l 5 >10 1 5 15 1 0 5 e 0 B ")( IS .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. 52 ....... )CEFVO\ I .. .. ...... 0 2 8 a 1 0 12 ..,,,,,._ ...,.,,.,. A .\h.:~ ... ..,.,, ...... Yi"RAOR V r'WII., ......... .. ~. . 0 2 8 a 1 0 12 .:,Z:.,M /\ 2EAAQA ...... .. ,., .... 0 2 4 6 8 1 0 12 Time (sec) Figure 3. EKF position errors with cr s tandard deviations. JX'YAA >l(VAA XEIVIOR , .. ......,........ ... ". 0 2 6 8 1 0 1 2 .,,.., .. 'rvAA VERROR \ . . . 0 a 8 10 1 2 ffiAA r.. ,,,,.,. \ __ ................. z:eAAQA .I . . 0 2 a e 10 12 Time (sec) Figur e 4. SDREF position errors with a s tandard deviation s. +3XVAR 3XVAR XERR O R .. " . : 1 I  0 2 a 8 10 12 .JWAR YVA R A AIJ' .... YERROR . ~ I 0 2 8 a 10 12 .:Jr./AR I\ .JZVAR ('\ l \ _.,...., __ .. .......... ___ ~.. ......... ZERROR  .. .. ,I . 0 2 a 8 10 12 Ti m e (sec) Figure 5 Boot s trap po s ition e rrors with crst andard deviations.
PAGE 64
15 10 5 E 0 Cl) 0 a.. >5 10 15 15 10 5 E C: 0 ;P CJ) 0 0 a.. >5 10 15 .. .. .. .. .. .. .. .. 1.0 5 3 EKF EKF .......... SDREF .......... SDREF ........ . . . . . ........ ....... ... . 0 2 4 6 8 10 1 2 Time (sec) F i g u re 6. Compar i so n of EKF and S D REF covar i ance cr sta n dar d deviat i o n s for Y pos iti o n 0 5 EKF E K F ............ SDREF ........... SDRE F .......................... .......................... .. ......................... ........ ..... . ....................... ......... .................... ~.... ~ I 0 0 0.5 Time (sec) 1 0 1 5 2.0 F i gu r e 7. Ex p anded v i ew of F i gure 6.
PAGE 65
54 6.2 PDF Distribution Comparison Using the baseline scenario, a comparison of the posterior (afte r the measurement update) probability density functions was investigated, to determine whether parameterizing versus l inearizing the nonlineariti es has an effect on the shape of the pdf. The scenario given in Table 1 was run once and the s tate estimate and stan dard deviation were computed three times during the engagement. Sampl es were taken at 1, 5, and 9 seconds. The pdf's for the EKF and SDREF were computed using a Gaussian distribution function based on the mean and stan dard deviation at each time increment. The resulting curves, Figures 8 through 9, do not represent the ''true'' pdf as nonlinear dynamics are involved. However they do s how the ba s i s of the theoretical filter f orrnulation. (''True'' pdf functions, calculated u s ing Monte Car l o sim ulation s, are discussed in following paragraph s.) The pdf for the bootstrap fi l ter was generated by comput in g a hi stogram of 1 ,000 samples, curve fitting a function u s ing a c ubi c sp line and nor1nalizing the area under the curve. Figures 8 throug h 9 show the oneru n posterior pdf functions for the filters. The first time increment is taken I second into the flight. The filters are still trying to converge and therefore the distributions are large. Thi s is demon s trated in F i gure 8 by comparing the narrow spike of the highly converged bootstrap filter to the wide, non converged pdfs for the other filters. If consideration is given to the scal ing on the xaxis, no concern shou l d be g i ven to the fact that the di stributions are not near the tru e s tate. It can be see n that th e estimates are close to th e ''truth''. The xpos ition state shows a more acc ur ate representation of the pdf than theyand zpos ition states. This is due, in part, to measurement noise decreasing as the engagement progre sses, and a more accurate
PAGE 66
........ N ..... X a. ........ "' N 55 True St a i 0.025 ... Bootstrap T ime 1 Second EKF .. .......... SORE 0.020 ... .. 0 015 .... .. 0 010 ._ .. 0.005 ._ .. 0 .000 .~ .. ~ :::.:.::::. . ... .. i,. iv. _ ........  . ,'' ... :',':.:; ...... ... .. "" ..... ,._ ..... :...,i~,,;,,. .... ......... I I I I I I 88500 89000 89500 90000 90500 91000 X P osition (m) 0.05 ., .. 0.04 ._ .. 0 .03 ... Time = 5 Seconds True ,tat e Bootstrap EKF .......... SDREF 0.02 .... c a. CJ) N X a. 0.01 0 .00 0.16 0.14 0.12 0.10 ... .. 0.08 ._ .. 0.06 ._ 0.0 4 ... .. 0.02 ... 0.00  I I I I I 50000 50200 50400 50600 50800 51000 X Position (m) Time = 9 Seconds ~ . . .. "" .. ~,; I, .... ... .. __ .. ...... ;., ......... True Sta . ~. . '. ..... .... . '. .... ... .. J ,_ Bootstrap EKF .......... SDREF ... .. ... ... . . .. . .. ... ........ ::.:0. 02 'l1'''L''__ ...__,~..._ 'iL'&........'..l''L'' 12180 12200 12220 12240 12260 12280 X Position (m) Figure 8. P osterior probability density functions at tim e = 1 5 and 9 seconds for X position.
PAGE 67
0.25 0 .20 0 15 N >,. 0 10 ......... a. 0 05 0 00 0.5 0.4 0 .3 .,, N 0 2 >,. ......... a. 0 1 0.0 5 4 3 0) N 1 0 True Stat T i m e = 1 Second / . 56 . . .. , ', ' ,. '. . . '. , , .. .., ,_.:. .. .... ....... ' . ' , ' , , . ' . \ . I ' \ ' ' ' ' Bootstrap EKF ..... SORE 27 4 272 270 268 266 264 262 260 258 256 Y Posi tion ( m ) ....... .. .. . t ,. I : ,, : \ Time = 5 S e con d s I : : .' : : : : : ,' :' . : l/ I .' ': : : I ': I,' I; ,: : .. .. .. ,. '. ,. True S tate \ ,1 ' I . \ . \ . \ I I \ I I , \ I , I , ' , I ', ', Bootstrap EKF .......... SDREF '""" ...... .. .. .. .. . ~421 420 419 418 417 Time= 9 S econds Y Posi t io n ( m ) I I I \ I \ .' , f .. I o t : : I 1 I ,; 1 I J :' I ' : : I I : l : : ' I I I ; \ ; : I l ' ' : J ' ' : I ' I ' ' ' . r . . . . . ' ' I ' ' ' ' : ' t . \ ' . ', 4 16 4 15 Bootstrap ....... EKF ,_, ...... SDREF True State ... ... . ... .0 173.5 1 73.0 172.5 Y P os it ion ( m ) Figur e 9 Pos t e rior prob a bility den s ity function s at time= 1 5 and 9 se conds for Y po s ition
PAGE 68
57 knowledge of the true state is obtained. Another interesting factor shown in Fi g ure 9, is t h e compression of the EKF and SDREF distributions. Thi s i s typic a l of the "fal lin g asleep" problem associated with many filters. As th e s tate estimation error become s s mall the covariance matrix al so be co me s sma ll, and the filter react s slowly to new data, wh i ch takes the estimate away from the converged answer. The boot s trap filter h owever, keeps the distribution open wider. Even u s in g the onerun pdf' s, the EKF and SD REF pdf s hapes appear to be similar. To determine the ''true'' pdf functions for the EKF and SDREF the filter s were run in a 100run Monte Carlo experiment Filter state estimates were save d at the same time increments. The s ample s from each of th e run s were u se d to generate th e pdf curves in Fi g ure s 10 through 11. The same hi stog ram method described in th e previou s paragraph for the bootstrap filter was u se d t o plot the curves. ( The boot s trap filter did not need to be examined in this fa s hion as it s forrr1ulation i s inherent l y Monte Carlo.) Th e largest differen ces can be see n in Figure 11. In the later time s ample s (5 and 9 seco nds), the previously tight distributions have flattened out, while the wide distribution s hav e tightened up Thi s is to be expected as the onerun case is only one realization of th e Monte Carlo runs. Summari z in g the above ana l ysis s how s that the EKF and SDREF pdf' s do n ot re se mb l e the ''tru e'' pdf. How eve r they are s imilar to each other Thi s may b e ex pl a ined by the fact that the EKF linearization method provides an accurate filtering methodology as long as the filter s tate e st imat e is close to th at of the ''true'' s tate. If a tar ge t maneuver is introduced that causes the state estimate to be far from the ''truth '', then it is a ss um e d
PAGE 69
N .... X 0. 0 .030 0 025 0 020 0.015 0 010 0 .005 0 000 58 True State Time = 1 Second . ; . ... . . ' : .. ... : ,, "''. : .' !, . ' : ,I \ ... \ / .' I / ' ' : \ I , \ ,' I \ : I \ , I , \ t : \ .. : ' : I ". \ : , .. \ .,, \ :, . J:, ... \ I , ,f \ '. ' ~\ I/ . . . . ' Bootstrap EKF ......... SORE ~# ', ... .. ... .. 89720 89740 89760 89780 89800 89820 89840 89860 X Position (m) 0 .05 , "' N X 0. 0, N :! X a. 0 04 0 .03 0 .02 0.01 0 00 50440 Time = 5 Seconds True late . . . . . ' ., ' \ : .. : \ .. \ . : , . \ : : \ r ', , I I : I I I : ' t I I : \ : . .:,' ~, \ .: ,' . I \ \ . \ / I , I .. t : . : ,I ' : .. \ : . :, . .. ,: 50460 . ...... 50480 50500 50520 X Position (m) Bootstrap ............ EKF .. SDREF 50540 0 16 . 0 14 0 12 0 10 0.08 0 06 0 .04 0 .02 0.00 Time = 9 Seconds ,.:, True Sta I .'t I _. ' t .. ; : : : ' : I : ' : I : : ' : : : . : : . ' .:, : : ' : .. , , ... \ ': I \,  Bootstrap  EKF .. SDAEF 0. 02 ,__ ____ __. _____ __._ _.____ ___ ....._ ____ 1 12220 12240 12260 X Po si tion (m) Figure 10. Monte Carlo posterior probability density functions at time= 1 5 and 9 seconds for X position.
PAGE 70
59 0.35 ....., N >, 0. U) N 0.30 0 25 0.20 0.15 0 10 0.05 0.00 0 5 0.4 , 0.3 ., ., 0 1 ., 0.0 ., True State . Time = 1 Second : ,, I . . . . . . . . . . ...... __ .... 270 I I I I I ' ' ' I f \ ' I , I \ ., ,: : I \ I I . ' . ' ' . . ' I I I I ' I ' I \ \ Y Position (m) T rue State I' Time = 5 Seconds .~ _,.r ' I _,.. .. r  ,,, .. .,,r ,, ,r r ... ~ I I I 425 420 415 Y Position (m) .., ~ ., Bootstrap EKF SORE 260 Bootstrap EKF SDREF . .... I 410 405 2."' N 1 0 Time = 9 Seconds 175.0 174 .0 Bootstrap EKF ..... .... SDREF True State 173.0 172.0 171. 0 170.0 Y Posit io n ( m ) Figure 11. Monte Carlo posterior probability density function s at time= 1, 5, and 9 seconds for Y position.
PAGE 71
60 that the distribution s between the EKF and SDREF differ dt1e to the different theorie s of linearization versus parameterization. Th e following paragraph s explore thi s sce nario. The filter s were once again run in a MonteCarlo s imulation However this time the target perform s a thru s t c utoff maneuver, 2 seco nd s a fter launch which cau ses the target acceleration in all axes, to drop to zero for the re s t of the flight Thi s is repre se nt a tive of timing the engagement too c lo se to booster burnout. The filter s are s till operating under th e assumption of a GaussMarkov target acceleration model even though no target acceleration i s present aft e r 2 seconds. The same initial conditions liste d in Table 1 were repeated Figure s 1 2 through 14 s how the divergence from the ''true'' acceleration s tat e for the filter s For previou s ly di scusse d rea so ns, Figure 14 s hows the boot s trap filter perfor1n s poorly when tracking the acceleration s tates. The deviation from the ''tru e'' s tate i s signif icant after the maneu ve r occurs Thi s i s a g ood test to deterrnine th e differences between linearization and parameterization Once again the Monte Carlo pdf functions were plotted in Figures 15 through 16 Sample s were taken at 2, 2.5, and 3 seconds into the engagement This allow s the pdf to be examined before and during the mo s t s tre ss ful tim e o f the man e u ver. Supri s in g ly, th e EKF and SDREF pdf s h apes d o not differ significantly. It can be see n in Figure15 that although the ''true '' probability ha s become bimodal, the s hapes are s till s imilar While the s inglerun EKF and SDREF di s tribution s can n o t be bimodal Figures 15 through 16 repre sent Monte Carlo res ults. The result s from the pdf co mpari so n experiments indicate that there is n o appreciable difference in pdf shapes between the EKF and SDREF The conclusion may
PAGE 72
0 ao l 60 <10 20 >0 eo <10 l 20 0 N <10 ..a .. .. .. .. .. .. .. I .. eo eo <10 20 0 eo ao <10 20 0 80 <10 20 0 20 o ,0 .. 1 .. .. .. 1 .. .. I ao 15 50 N 2S 0 N ~ .. .. .. . 61 ____ """" ""' . . . 0 2 4 e 10 T TAAM "' ', '\ ',. ./ ... 0 2 4 e 10 __ _,...,. ....,/ .r,.._ .,/ . 0 2 e a 10 Time (sec) Fi g ur e 12. EKF es t i m a t e d a nd tru e acc el e r a tion SXT AXREAL 1 .....'.. ,.,r ........ ,,_.,,,.... _ ~~. \'"\...._ i.... ... ~. . .""lo.. 0 1 2 3 5 7 a 10 11 SYT 1\,"'""'.1""".,, r~ AYREAL \ .J ,. . . . . 0 1 2 3 5 7 9 10 11 SZT AZREAL ._ .. . .,.., .......... . . . 0 I 2 3 5 e 7 9 10 11 Tlme (sec) Figur e 1 3. SDREF es tima ted a nd tru e a cce l e r a ti o n ... .. . ............ . .. .. ... .. ... . 0 2 i .. .. . .. ... . 0 2 .. .. . 0 2 . XN '"' .. .. ... . . a a 10 = t 1,... _.,._ ... __ A..i. . .. .. .. .. .. . .. . e 8 10 .... z:n, .. . ... ,r ... ... ' 5 Tlme (sec) . 8 10 12 12 12 F i gu r e 14 B oo t s t rap e s ti m at e d a nd tru e acce l e r a tion
PAGE 73
""1 N ,:: X a. 0 .030 0 025 0 020 0 .015 0 .0 10 0 .0 05 0 .000 0.02 0 .01 .,; N Ill X ... Tlme = 2 Seconds ,. ... ., l ,. ,: " h ,: .: ,. ,; J' ,. .. 62 \ I < I \ ', \ \ \ \ \ \\ \ ,, B oo tstrap EKF ....... SDREF , ".__ . 79820 79840 79860 79880 79900 79920 79940 X Position ( m) Time 2 5 Seconds Bo o tstrap T rue State .. ... EKF ~ SDREF .......... . . . .. .. .. :, . ( : : . \ I : i . I l : t ' ' : .. : ' :, : : ; .. : . :, I h I r t i I I > l' I l : ' / ,. : l . en N X a. a. I : I : ; : ' 0.00 ~ .r \ ..... ,. .. .. I 74800 749 0 0 0 .020 T ime = 3 Seconds r 0.0 15 ,. .. .. .. . . : /l .. > I ' \ ., . 0 .010 .. I ,; ,: ll ., I, \ ,: ': ,, 'I ., l '/ : ' '' ,, 0 .005 .. ,: ,.  : \ .. : i ,: ,. '' . .. ,. : ,. I I 75000 7 5100 X Position (m) T rue Sta : . .. : I I \ I : : I ~ I 75200 75300 B oo t s trap EKF . . . iA~ :, 1 : .. .. f \ I \ / I ..... S DREF \ / i I ? I .. I i I ; .. \ .. 0 .0 00 ,... ~\ \ .,!r ,., ~ ~__.c::,__, .. ..., I I I I 69900 70000 70100 70200 70300 X Posi t ion (m) Figure 15. Monte Carlo posterior probability density functions at time= 2, 2.5, and 3 seco nd s for thru s t cutoff X posi tion.
PAGE 74
"' N >. ......, C. "' c,j N "1 >. ......, C. (') N (') ,_ >. ......, C. 63 0.35 .,, 0.30 0.25 0 .20 0 15 0 10 0.05 True late Time = 2 Seconds .. .. ' I . ' ' / ' .. . . .... . Bootstrap EKF SDREF \ '1 . . . ' . 0.00 390 380 0.4 0.3 0.2 0 1 0.0 Time = 2 .5 Seconds Y Position ( m) ... .. l j .. ' .. \ . : : ' l \ ' : I : True S tat e . ' ' . 405 . ' I I I ' . ' . . , ' I \ ' ' ' ' ' ', 400 Y Position (m) Bootstrap EKF SDREF 0.40 r,, 0.35 0.30 0.25 0.20 0.15 0.10 0 .05 0.00 404 Time = 3 Seconds ... ,. .. ... 402 ., . \ ~, . \ .. , ' \ :, ~ .. I .. .. ' # : ' : : ' : : I : ' : : : I : I : : I / 400 398 Y Pos ition (m) Bootstrap EKF SDREF True State 396 39 4 Figure 16. Monte Carlo posterior probability density functions at time= 2 2.5, and 3 seconds for thrust cutoff Y position
PAGE 75
64 be drawn that the two filter s operate with identical statistical models although the filters' method of handling and operating on the stochastic processes are significantly different 6.3 Sensitivity Analysis The following section examines the SDREF sensitivity to measurement noise, initialization error and the robustness properties in the presence of stressing target maneuvers. Only the EKF and SDREF are used for this area of the study. Further bootstrap filter development is needed to make it useful in less than ideal conditions The EKF and SDREF were tuned separately to give them the best chance of performing w e ll under these varying conditions. Therefore each filter may have been started with different initial conditions and statistical information. The following studies are accomplished using perfect state information to guide the vehicle, while the EKF and SDREF run open loop. Allowing the guidance algorithms access to perfect state information removes the closedloop guidance algorithm effects from the analysis of filter performance 6.3.1 Initialization Error To test the response to initialization errors, the filters were given erroneous initial state estimate error s of 1 and 5 percent The EKF initial covariance was set to approximate this e1Tor. The SDREF, however has no direct means of initializing the covariance matrix as it is solving for the covariance matrix at each filter update time. This will be shown to be a detriment for this particular implementation of the SDREF. Table 2 gives the initial filter conditions for the EKF and 1 percent initial error Measurement noise was held constant for both cases, as was the measurement noise matrix given to both filters. Range measurement noise is again range dependent and i s
PAGE 76
65 given as a percentage of ''true'' range The initial setup for the SDREF used the same parameters as those of th e EKF except that higher values of process noise were tried to increase the covai iance respon se. The only tuning parameter s in the SDREF are Qk_1 R and the nonlinear parameterization methodology The measurement noise matrix R is not directly related to the i nitialization error; therefore only the other two tunin g parameters can be used. Because this implementation of the SD REF bas linear dynamic equations the only method available for tuning i s to adju s t Q k i This is an unacceptable way of tuning the filter for initialization e rrors as it causes the covariance estimate to remain artificially high fo1 the entire flight. The SDREF as previously derived doe s not permit any initialization of the covariance mat1ix a s it i s simply calculated from the s olution of the SDRE. Thi s deficiency is examined later in thi s sec tion Tabl e 2. EKF and SDREF parameters onepercent initialization error Parameter s ( units) X y z Aspect Angle ( degree s) 30 30 30 SIGP ( m ) 1000 35 30 SIGVP ( m l s) 100 3 4 2.2 a 2 r ( m 2 ) l .OE6 l .OE3 1 .0E3 SIGXHAT a 2 V ( m 2 / s2) l .OE4 50.0 50.0 a 2 A ( m 2 I s4) l OEl 1.0El 1.0El R ( m ) l .OE5 10 0 10.0 XHAT V ( mis) l .OE4 335 9 223.9 A (ml s 2 ) 55.2 57 4 38.3 SIGQ (EKF) (m2 / s4) 196.0 196.0 196 0 SIGQ (SDREF) ( m 2 / s4) 1.0E6 1.0E3 1 0E3 Meas Noi se (%,m,m) 0 1 50.0 50.0
PAGE 77
66 Figures 17 and 18 show the positionstate response of the EKF and SDREF to a Ipercent initial error. (Note the oscillations in the SDREF estimate during the initial settling portion of the flight.) Velocity and acceleration trends were similar. The filters were then run with a 5percent initial error. (The initialization parameters for the EKF and SDREF are given in Table 3.) Because of the type of engagement and large distances and closing velocities involved, a fixedpercentage type of initial error induces large disturbances in the xdirection. Therefore, little attention is given to the performance in this direction, as these magnitudes of error are not likely to occur Figures 19 through 20 show the performance for a 5percent initial error. The inability to properly adjust the initial covariance matrix in the SDREF causes a large oscillatory error in the early portion of the estimation problem. This error was unacceptable and a new initialization procedure was investigated The basic problem with the initialization of the SD REF is that the initial covariance matrix is computed using the solution of the algebraic Ricatti equation. This indicates that no large errors are expected. Therefore, there is no method of supplying the SDREF with an initial covariance error of the magnitude of the expected initialization errors. Because of this the SDREF computes a smaller magnitude covariance matrix and therefore a smaller gain to be used in the state update This is why the filter oscillates during the initial portion of the flight. If the SDREF could be "jumpstarted" with a high initial covariance, the majority of the initialization error could be taken out in the first time step.
PAGE 78
300 .. 200 .. 100 E 0 g Q. 100 X .. .. .. 200 .. 300 10 lo 0 lo 10 10 lo 0 1000 .. 800 lo 600 ... 4 00 200 X 0 lo 200 400 40 .. 30 .. 20 E g 10 Q. .. >.. 0 .. 40 lo 30 20 E 10 g Q. .. .. N 0 .. .. 20 0 0 0 p ..... .,,.. .., ........ . . . / 2 4 \. :: ... to., .... :" .oi .. "' ......... ~ ... ..... .......... I 2 4 \. ,. .... ... ' ., . ........ ... ,, ... .......... ,, .......... / 2 4 ; 67 ...... 6 ~.. .. 6 .. ... 6 Time (sec) P1<3PLUS PKJMINUS XHATEAA 8 10 Pl<3PLUS PKJMINUS XHATEAA .. ........ .. 8 10 PK3PlUS Pl<3MINUS XHATEAA ... .. 8 10 Figure 17. EKF position errors with astandard deviations for I percent initialization error. +3XVAR ..SXVAA XEAROR ~, ,, ,. V 0 2 4 6 8 10 +3VVAA VV AA YERAOA . J ~.,.. 0 2 4 6 8 10 +3ZVAR ZVAR ZEAROR ,,... ..... r,,..,.,...._ \ / ._ ,_ .J . 0 2 4 6 8 10 Time (sec) Figure 18. SDREF position errors with astandard deviations for Ipercent initialization error. 12 12 12
PAGE 79
400 900 200 E 100 0 c.. X 100 200 400 10 I" I" 10 10 .. 0 I" 4 000 .. 3000 .. 2000 E .. c.. 1000 X 0 .. 1000 1 50 .. 100 e 50 .. 8 0 c.. > .50 I" 100 I" 150 100 80 .. 60 I" 40 .. e 20 l 0 N 40 I" .. .. .. .. I" I" 0 . : " ..... ~.. \. ., .. : ':. .:~, ... '. 2 .. .. \ . ..... , 0 2 ,. t ... , :' ~ ~ i.. .. ..: ....... r 0 2 4 .. ....... _, ... 4 .. ,.; .... 4 68 .. 6 .... 6 I., 6 Time (sec) .... PK3PlUS PK3MINUS XHATERR 8 10 PK3PLUS PKSMINUS . XHATERR .. 8 10 PK3PLUS PK3MINUS XHATERR .. 8 10 Figure 19. EKF position errors with astandard deviation s for 5percent initialization error. .:JXVAR XVAR  XERROR .J9.... \ /' ; .. 0 2 4 6 8 10 +3VVAA 3YVAR .......... _., ERROR \ \ ,,, .... " . 0 2 4 6 8 1 0 +3ZVAR 3ZVAR ...... .. ZERROR .. / ,_.,. . 0 2 4 6 8 1 0 Time (sec) Figure 20. SDREF position errors with astandard deviations for 5percent initialization error. 12 .... 12 ,.. 12
PAGE 80
69 Table 3. EKF and SDREF parameters fivepercent initja]ization error Parameters ( units ) X y z Aspect Angle ( degrees ) 30 30 30 SIGP (m) 5 .0E3 168 0 112.0 SIGVP (ml s) 5.0E2 16.8 11.2 (jr 2 ( m2) 25.0E6 3.0E4 3.0E4 SIGXHAT (j2 ( m 2 / s2) V 25.0E4 3.0E2 3.0E2 (j 2 A ( m 2 /s4) 1.0El 1.0El l.OEl R (m) l.OE5 10.0 10.0 XHAT V (m/ s) 1.0E4 335.9 223.9 A (ml s 2 ) 55.2 57.4 38 2 SIGQ (EKF) (m2 / s4) 196.0 196.0 196.0 SIGQ (SDREF) (m2 / s4) 1.0E6 l.OE3 l .OE3 Meas Noise (%,m,m) 0.1 50. 0 50.0 If the gain calculation used in both the EKF and SDREF is examined, it can be seen that the covariance used at time step zero, ( ) is needed in the calculation. If it is assumed that the value of P k that the SDREF computes is P k ( +) (i.e., after the updated measurement has occurred), then the SDREF is initialized with P k () to be used at time s tep zero. Using this methodology the SDREF was initialized with the same covariance as the EKF and the SDREF process noise was also reduced to the same level as that of the EKF. This change in methodology had positive impacts on the performance of the SDREF when confronted with initialization errors. Figures 21 through 26 show the new perfor1nance of the SDREF using the same initialization errors as in the baseline runs. It can be seen that the initial oscillatory motion in the state estimate errors has been corrected and the filter compares favorably with that of the EKF.
PAGE 81
70 100 ..,..,.,. .. .:IXVAA ,,,.~ ..'\, lCEAA0A 0 ~ ,/ . ~ '" X I t .. > )( i CD > > I ai > N 0 2 a 8 10 5 .. 'IYVAA ........ YM0f\ l 0 V\,.,1" '\ .I\/', I\ N\ A,. .. ~ vv~ v _....,..,,,_ > 0 2 8 a 10 e s 4ZVAA 4 '/NAA 0 l'A ,. ,.A, ., ... .... "\. A . .. N i )( "' i > i N .. 0 4 8 10 Time (sec) Figure 21. SDREF pos ition errors with astandard deviation s after new initialization method, Ipercent error. 200 'SI/WAA 100 _/v,..__..._ 'N'l:'(AA ' 0 .. ''~ ~,..., .... ~ .. .. .. 100 20 0 2 4 8 8 10 15 10 5 0 .5 10 4'('(\fA/1 .... ~VYERR A i' '\ ~ _,, A a \1\1'1 ..... 'I.. A ,. .... .&I'\ ; \. ' , L V .. 15 20 0 2 4 20 8 8 10 0 I 0 2 8 a 10 Tim e (sec) Figure 22. SDREF velocity erro r s with astandard deviation s after new initialization method Ipercent error. &O 4() """""""  20 .. A Al(EM 0  ~,, ' ...... ./ .... ....... I" '\. ., "" 60 0 2 4 40 6 8 10 .. .. :IM'\IAR 20 ""''"'""  .. AVEAPI .. 0 ,., r V . "'"'I. ; ~ ,,,,..\ti", '.I' 'vN" . ....., .. r .. 40 . 4() 0 2 8 8 10 .. ... 1'1 .... 20 0 ..... .. AZE~ .. I"',_ .J .I"'\ A.,. .J' .,, ~~.A;,,,. ... I\ I 1,./"" V . .. 20 I .. 40 0 2 8 8 10 Time (sec) Figure 23. SDREF acceleration errors with astandard deviations after new initialization method, Ipercent error.
PAGE 82
200 E IQO if 0 X 100 200 100 15 10 .. 5 0 > .5 10 15 15 10 .. 6 0 N 10 15 .. .. .. .. .. .. .. .. , 0 f / 0 / I 0 ~" ... 2 ' 2 .,. . 2 71 .m(VAR 4XIIA11 )(Mof\ .  e 8 to ....,..,. ./SYVAII vtiWIOR . . e 8 10 ..,.,.,. >ZiAA ZE 5 a 10 Time (sec) Figure 24. SDREF position errors with astandard deviations after new initialization method, 5 percent error. 1500 400 ai > 200 X 0 200 60 I 40 20 ai 0 > > 20 .40 80 80 0 40 $ 20 0 N 20 40 80 .. .. .. .. .. .. .. .. .. .. .. .. ..... ...... 0 2 ... . 0 2 0 2 ,_.,. 4V'I.VAA YXEAA . a 8 10 YVYVAA ~All  WEAR .... .... .. . . e 8 10 4'/l'tAA .:,vr,tAA VttAA ~ . 8 8 10 Time (sec) Figure 25. SDREF velocity e rrors with astandard deviation s after new initialization method, 5 percent error. .;; 50 0 i X 100 150 N 20 0 >40 20 !l 0 N 20 .. .. .. .. 0 .. .. .. .. 0 .. .. I .. .. 0 ..... . _..,.~'" 2 "'\_ ..., ', ../' 2 4 ..... ......,___ . ., 2 ,11/,J!NAA 41>XVAA ............ .. All.,.. .. 6 a 10 'YVM .JAVVAR .. AYAA .. l",'.r....,,_,. '~ .. "..__~ ..... e 8 10 .UYAA 0JANAA AZEARM ... .,. ... ""',"~ ,... ,_~ .... ,,I e 8 10 Time (sec) Figure 26. SDREF acceleration errors with astandard devi a ti o n s after new initialization method 5 percent error.
PAGE 83
72 This new initialization methodology not only improves the perforrnance in the presence of initial errors, but also allows a convenient way of initializing the SDREF. 6.3.2 Measurement Errors This section investigates the performance of the EKF and SD REF when subjected to various degrees of measurement noise. The relative perforrnance of the two filters using the baseline measurement noise was shown in Figures 3 through 5. Using the same baseline initial conditions given in Table 1, two other levels of measurement noise were added (Tables 4 and 5). Table 4 EKF and SDREF measurement noise parameters level I Parameters Measurement Noise (%, ra
PAGE 84
73 The second level of noise is large for a typical engagement; however it exercises the robustness of each filter. Both filters were given measurement noise covariance matrices that matched the noise characteristics given to the sensor. The position state estimate plots are given in Figures 27 through 30. There is not a significant difference in filter performance using the large noises. As expected, the estimates for both filters are noisier and more erratic. The acceleration state estimate has become even worse under these conditions. This estimate was initially poor, and with noisy measurements the filters do not have enough inf orrnation to estimate the target acceleration to the desired degree of accuracy It can also be noted that the SDREF estimation error under the second level of noise seems to be worse than the EKF estimate. While these figures only represent one run, the SDREF seems capable of adequate perfo1mance even under stressing noise conditions. 6 .3.3 Target Maneuvers This section investigates the openloop perfor1nance of the EKF and SDREF against two stressing target maneuvers. The capabilities of the filters against a thrust cutoff maneuver have been previously discussed. In this section two other nonlinear maneuvers were tested. The first is a 10degrees per second rotating target maneuver. This maneuver allows the target acceleration vector to rotate approximately 110 degrees during the engagement. This is a moderately high target maneuver, but it is possible for an ICBM to perform this type of action in an active counterrneasure scenario. The second stressing engagement consists of a dogleg maneuver in which the target acceleration vector rotates 90 degrees out of plane and then reverses itself to rotate back inward 90
PAGE 85
E X E s Q. >E s Q. N E s Q. X g s 300 200 100 0 200 300 30 20 10 0 10 20 30 30 20 10 0 20 30 400 200 0 200 10 Q. 0 >10 10 5 0 .5 10 ,. ,. ,. ,. ,. ,. ,. ,. ,. ,. ,. "" "' "' "' ,. ,. "' .. ,. ,. ,. .. .... ... .. 0 2 4 '" :.. ... .. I' ~ I"' ... ....... . .. ... r 0 2 4 '" ,._ ...... ..... ,. .. .. ...... r~O' . .................... 0 2 4 . 74 .. .. .. .. . .. ... . . 6 .. .. .. 6 .......... ~ 6 T i me (s ec) .. .. .. . .. . .. . ~ . PKSPLUS PICSM INUS XHATERR s ..... . . .. ....... . 8 10 PK3PLUS PK3MINUS XHATERR . 8 10 PK3PLUS PKaMINUS XHATERR . 8 10 Figure 27. EKF position errors with astandard deviation s Ipercent error in X, 100 microradian error in Y and Z. 0 0 0 V\"~..1,'\ /',.,.. .... ../ ....... '.,.l'~ ........ r' \..\v,""\ 2 4 6 8 \ Ari\. .,t\ A. I"'. .. ..A ,, v......., \\ .. A ,.,1v ""\~ '" ,L ... .. "'\,I\ It" . ; 2 4 6 8 \r\., r,.I .. .I . 1/'""' ... t".. \ vr V v..~ ~I Vy V ,.._.,....,., . 2 4 8 Time (sec) 8 lXVAR XVAA )(ERROR . 10 .:JYVAR YVAR ......... .... YERAOR 10 +3ZVAR ZVAR ZERROR 10 Figure 28. SDREF position errors w i th astandard deviat i ons Ipercent error in X, I 00 microradia n error in Y and Z. 12 12 12 12 1 2 1 2
PAGE 86
1000 800 "' 600 "" 4 00 "' 200 E l 0 200 X 400 "' I" "' "' 000 "' 800 I" 150 "' 100 "' 60 ,. 8 0 0. .. >50 1 00 1 50 160 "' 100 ,. 50 E .. 8 0 0. "' N 50 "' 1 00 .. 150 2500 2000 1500 ,. 1000 I" E 500 8 0 0. .500 X 1000 I" "' "' 1500 .. "' 2500 "' 30 I" 20 "' 10 E "' 8 0 0. "' > "' "' 30 "' 20 "' 10 E "' fl 0 0. "' N 1 0 "' "' 30 0 0 2 \ .: '. .. ......... ... :f. r 0 2 \ ..... .... ..... ' I 0 2 75 .......... .. . . . .. M 4 .,_ .... .. .... . .. 4 ........... .. 4 6 6 6 Time (sec) PKSPLUS PK3M I NUS XHATERR ... ... ...... ........... .. __,,,r 8 1 0 PK3PLUS PK3MINUS XHATERR .. ........ 8 1 0 PK3PLUS PK3MINUS XHATERR . . 8 10 Figure 29. EKF position errors with cr standard deviations 10percent error in X, 500 microradian error in Y and Z. ..JXVAR 3XVAA XERROR \,J\ ...... .. ,,.. ......... __ ,. ,..,. ____ ..,,,,.. .. '.,,,....,., ...... ___ >< ... .,....... ,.,. ... . 2 4 6 8 10 +3VVAR Jo 3VVAA ;v\  VERROR f\ ., v 'Al \. vv ~ ..... ,. A 0 2 4 ,, .. ,.. v"'I "' ,n, T \. A ,,....,... 0 2 4 .I .... ""'I \,. \.I ""\ / 6 8 ,.~ ,.,...A'\ ,~\.J ., 6 Time ( sec ) 8 ~... y,. .. 1 0 +l!ZVAR 3ZVA R ZERROR ~;~ ... 1 0 Figure 30. SDREF position errors with astandard deviation s 10percent error in X 500 micro radian error in Y and Z 1 2 12 1 2 12 12 12
PAGE 87
76 degrees. The initial conditions used for this analysis were the same as the baseline scenario except for the process noise characteristics given to the filters (see Table 6). Table 6 Maneuvering target parameter s Parameter s ( units) X y z A s pect Angle ( degrees ) 30 30 30 SIGP ( m ) 0 0 0 SIGVP (mis) 0 0 0 a 2 r (m2) 1.0ES 1.0ES l.OES SIGXHAT a 2 V (m2 I s2 ) 1.0E3 1.0E3 l .OE3 a 2 A ( m2 I s4 ) 1.0E2 l.OE2 l.OE2 R (m) 1.0E5 10.0 10 0 XHAT V (mis) 1.0E4 335.9 223.9 A (ml s2 ) 55.2 57.4 38. 2 SIGQ (EKF) (m2 / s4) 866.0 866.0 866.0 SIGQ (SDREF) (m2 / s4) 1.0E3 l.OE3 1E3.0 Meas Noise (%,m,m) 0.1 50.0 50.0 The process noise levels were tuned to allow the filters to open up and "accept" the u nexpected target maneuvers While this allows for better tracking, it does lead to a noisier state estimate in both fi l ters. Figures 31 through 34 show the relative perfor1nance of each filter against these two scenarios. As can be see n from these figure s, there are no appreciable differences in the s tate estimates between the two filters They each track this maneuverable target, however the state estimate error is noisier and larger as expected. It should be noted that once the initial transient terms in the EKF covariance have sett l ed, they are similar in shape and value to the "steadystate" covariance terms of the SDREF.
PAGE 88
200 150 100 50 .. 8 0 a. 50 X 100 .. 160 to 1 0 ,. 0 ,. 10 1 0 E ,. Cl) 0 a. 0 N 80 ,. 80 ,. 40 ,. 20 E ., 0 0 a. 20 X ,. ,. ,. 10 8 ,. 6 4 E 8 2 a. 0 > ,. 4 ,. 8 6 ,. 4 e 2 8 0. 0 N 77 PKlPLUS PKlMI NUS \ lCHATERR '.. . ... . .. ,.... , ., .... ,. .. ..... .. . . ...... . ..... / I . 0 2 4 6 B 1 0 0 0 t. 1 ;, ..... :. ,.,. ""' ',.;~ ::t# "' .. .....  ~. .,.. ... .. .. ( ..J...,,;i,ii 2 4 .. ..... : .. . .. .. . .. ~.. ... .. ........ vi ....... 2 4 .. .. .. ... 6 6 Time (sec) .. PK3PLUS PK3M INUS XHATERR ...  8 10 PKSPLUS PKSM I NUS XHATERR 8 1 0 Figure 31. EKF position errors with astandard deviations for 10 deg/sec rotating target maneuver. +3X\IAR X\IAR XERROR Jv'\. "'\ ,. N..,, ,.,~..,I",. _,, ,_~ _... 40 .r ... ,.i .,/ V . .,,, .. .... "'v'V "J I" 0 2 4 6 B 10 +3VVAR 3VVAR YERROR _,A l''. I ~" . J ,ilA "' /ti l A.~ I \ ,, L. ft ., :i:i "\J' '\.1 I\ rV' '\~J' .. V "\.JI' y e..~, . 0 2 4 6 8 10 +lZVAR ZVAR ............... ZERROR ,.I"' .._,., I''\. _r .... .,..~ lA.. Al' ,. ,I\ ..,. v \ ~, . v, .. ~ "\ l\,,J 0 2 4 6 8 10 Time (sec) Figure 32. SDREF position errors with astandard deviations for 10 deg/sec rotating target maneuver 12 12 12
PAGE 89
200 .. 150 .. 100 .. E 50 .. 0 )( .. .. 100 .. 50 .. 200 10 .. 0 .. 1 0 .. 0 .. 1 0 80 .. 60 .. 40 .. e 20 .. 0 0 X .. lo 40 .. 60 .. 80 10 8 6 .. E 4 .. 2 0 >.. .. 2 .. 4 lo 6 8 .. 6 .. 4 e 2 8 c.. 0 N .. 2 .. 4 .. 6 \ '....... ..... ,_, ....... . .... ..... ............ ~ / I 0 2 4 ft ". ~,: :. .. h:4.4 j :\ ........ . .~ i ., '\ ft ... .. .. 0 2 4 : .,; . ' . ... ... : ... ~It,., .~./ \ IC 0 2 4 78 ... .........  . 6 ft ~ ', 6 ' .. 6 Time (sec) ... PK!IPLUS PK!l M I NUS XHATERR . 8 10 PICJPLUS PKSM INUS XHATERR ,.T . 8 1 0 PK3PLUS Pl<3MINUS .. XHATERR .. . 8 10 Figu re 33 EKF p os iti o n e rr o r s with 3 a s tand a rd devia t i o n s for d og l eg t arge t man e u ve r +IIXVAR SXVAR ... ... J
PAGE 90
79 6.4 SDREF Closed Loop Performance 6.4.1 Other Engagement Aspect Angles Up to this point the EKF and SDREF performance in an openloop environment have been examined. The following analysis examines the closedloop perfonr1ance of each filter against various engagement scenarios. This includes a variety of engagement aspect angles, target maneuvers, initial condition errors, and measurement noise levels In this section, however, performance in terms of probability of hit and percentage of fuel used to achieve the hit are measured. Measuring the percentage of fuel used, allows a comparison of how noisy the state estimates are, as increased error in the beginning of the engagement result s in fuel being wasted performing unnecessary maneuvers. This analysis was performed usin g 31 Monte Carlo run statistics. The first analysis focused on three engagement scena rios Up to this point a 30,30,30degree closing aspect angle was used. A headon, 90degree beam shot, and a tailchase engagement are now added. Figure 35 shows the three engagements and the cotTesponding simulation parameters are given in Tables 7 through 9.
PAGE 91
80 z (0,0,0) Aspect Angle Head On _______ x y z (0,0,90) Aspect Angle Beam Shot y z (0, 150 0) Aspect Angle Tail Chase Figure 35. Aspect angle description.
PAGE 92
81 Table 7. Headon simulation parameters Parameters (units) X y z Aspect Angle ( degrees ) 30 30 30 SIGP (m) 0 0 0 SIGVP (ml s) 0 0 0 02 r ( m2) l.OE5 1.0E5 1 0E5 SIGXHAT a 2 V ( m 2 / s2) l.OE3 1.0E3 l .OE3 (j A 2 ( m2 I s 4 ) 1.0E2 1.0E2 l .OE2 R (m) 1.0E5 10.0 10.0 XHAT V (ml s) 10000.0 335.9 223.9 A (ml s 2 ) 55.2 57.4 38.2 SIGQ (EKF ) (m2 I s4) 866 0 866.0 866 0 SIGQ (SDREF) (m2 I s4 ) l.OE3 1.0E3 l .OE3 Meas Noise (%,m,m) 0.1 50 0 50.0 Table 8 90degree beam shot simulation parameters Parameters (units) X y z Aspect Angle (degrees) 0 0 90 SIGP (m) 0 0 0 SIGVP (ml s) 0 0 0 a 2 (m2) l.OE5 l.OE5 l.OE5 r SIGXHAT a 2 (m2 I s 2 ) l.OE3 l.OE3 1.0E3 V a 2 A (m2 I s4) l.OE2 l.OE2 1.0E2 R (m) l .OE5 10. 0 10.0 XHAT V (mis) 1.0E4 498.2 0 0 A (ml s2 ) 0.0 88.3 0.0 SIGQ ( EKF ) (m2 / s4) 196.0 196.0 196.0 SIGQ ( SDREF) ( m 2 / s4) 196.0 196.0 196 0 Meas Noise (%,m,m) 0.1 50.0 50.0
PAGE 93
82 Table 9. Tailchase simu l ation parameters Parameter s ( units) X y z Aspect Angle ( degrees ) 0 150 0 SIGP (m) 0 0 0 SIGVP (m/s) 0 0 0 (]"2 ( m2) l.OE5 l.OE5 1.0E5 r SIGXHAT (]"2 (m2 / s2) l.OE3 l.OE3 1.0E3 V a i A (m2 I s4 ) 1.0E2 1.0E2 l.OE2 R (m) 1.0E5 10. 0 10.0 XHAT V (ml s) l.OE4 0.0 237.9 A (m/ s 2 ) 76.5 0.0 44.1 SIGQ (EKF) (m2 / s4) 196.0 196.0 196.0 SIGQ (SDREF) ( m 2/s4) 196.0 196.0 196.0 Meas Noise (%,m,m) 0.1 50.0 50.0 The 90degree beam shot i s the most stressing of these engagements as the entire target accele1 ation profile is out of plane. The headon case is one of the easier engagements when using a range measurement, and the tail chase lies inbetween these two in terms of estimation difficulty. Table 10 displays the results from the Monte Carlo analysis against the se engagements. Probability of hit (PHIT) is the probability that the interceptor hit the target. Circular error probable (CEP) is the probability that 50percent of the misses lie within this bound 95percent mis s includes ninetyfive percent of the misses, and worst miss displays the worst of the Monte Carlo runs. The mas s usage numbers reflect the same quantitie s except in terrns of percent of divert fuel used. The vehicle has a mass of 4 kilograms, 1 5 of which is devoted to divert fuel. When this limit is exceeded, the interceptor will continue to fly towards the target with no further guidance corrections, typically re s ulting in large mis s distances. Since the se engagements
PAGE 94
83 had good initialstate estimates and the target is not maneuvering, all of the run s were hits; however as shown in the following sections that this is not always the case. Table 10 Miss distance performance for varying scenarios o o 90 150 Parameter ( units ) EKF SDREF EKF SDREF EKF SDREF PHIT 1.0 1.0 1.0 1.0 1.0 1.0 CEP (m) 0.185 0.262 0.432 0 .62 5 0.315 0.265 Mean Miss (m) 0.213 0.272 0.546 0.667 0.347 0.354 95 % Miss (m) 0 364 0.426 1 .2 12 1.209 0.654 0 .67 1 Worst Miss (m) 0.421 0.599 1 .828 1.679 0.803 0.895 Mass CEP (%) 16.70 27.61 67 .38 72.67 43.42 48 .06 Mass Mean (%) 16 90 28.27 66.98 72.76 43.11 48.40 Mass 95% (%) 21.03 32.78 69.25 75.22 46.10 51.62 Mass Worst (%) 22.70 33 16 71.20 75.51 48.57 52.59 From Table 10 no significant difference in the closedloop perfor1nance of these two filters under the above scenarios i s s hown All performed adequately and ensured a hit against the target in all cases. Fuel usage was also not a deci s ive factor as both filters differed by only a few percentage points. 6.4 2 Stressing Target Maneuvers To further test filter perf 01mance, two other target maneuvers were considered. These maneuvers are on the outer edge of filter capability of properly tracking the '' true'' state. They were included to inve s tigate whether the parameterization of the SDREF allows an advantage over the linearization of the EKF in highly nonlinear environments. The first maneuver is a rotating target however in thi s engagement the engagement aspect angle was a 90degree beam shot and the rotation rate was increased to 20 deg r ees
PAGE 95
84 per second Therefore the target will have completed an approximately 200degree turn during the length of the engagement The second maneuver consisted of a dogleg evasive maneuver. This maneuver is similar to the previously described dogleg maneuver except that it was performed in the 90degreebeams hot scenario, causing it to become more difficult for the filter to estimate Table 11 displays the results of Monte Carlo runs The process noise levels f or each filter were tuned to give the best performance against these scenarios Table 11. Miss distance performance for stressing target maneuvers Parameters (units) 20 /Sec Dogleg EKF SDREF EKF SDREF PHIT 0 935 0.806 0.871 0 742 CEP ( m ) 1.61 2.20 1.14 2 .85 Mean Mi ss ( m ) 1 .83 200.20 1.50 3.05 95% Miss ( m ) 4 30 7 44 3 94 6.95 Worst Miss ( m ) 6.00 6121.89 5 .59 7.36 Mass CEP ( % ) 100 100 100 100 Mas s Mean (%) 99 94 99 .99 99.87 99.98 Mass 95 % (%) 100 100 100 100 Mass Wo1st (%) 100 100 100 As shown in Table 11, the SDREF perfo1r1ed slightly worse against the targets than the EKF The one large miss distance for the SDREF wa s caused by poor tracking, which caused the interceptor to run out of fuel early and coa s t to a large mis s distance The purpo s e of this section was to s tress the filters to their extreme operation a l limits and inve stig ate the behavior. These type s of target maneuvers would be con s id e red
PAGE 96
85 unlikely. Howeve r they do provide a convenient way to detennine the robu s tness of each filter design Even at these extremes, the SDREF was almost able to meet the EKF performance meas ures 6.4. 3 Me a s ur e m ent E rr ors This paragraph discusses filter sensitivity to various levels of measurement error. For this analysis the range measurement error was held constant at 1 percent due to the s ize of the range error when dealing with large engagement distances. Azimuth and elevation noise s were included for 100 500, and 1000 microradian e1rors. The results displayed in Table 12 were done u s ing a 31run Monte Carlo experiment. Table 12 Miss di s tance performance for v arying measurement noi se Parameters (units) lOOrad EKF SDREF 500rad EKF SDREF lOOOrad EKF SDREF PHl 'I' 1.0 1.0 0.968 1 0 0 935 0.806 CEP ( m ) 0 455 0 .460 1.39 1 41 3.09 2.22 Mean Miss (m) 0 468 0.568 1 47 1.45 3.17 3.14 95% Miss (m) 0 .775 1.218 2.47 2 .60 5.66 6.29 Wors t Miss (m) 0.935 1.569 2.55 3.08 6.03 8.18 Mass CEP (o/o) 77.13 76.86 84.19 81.83 89.79 86.99 Mass Mean (%) 77.57 76.96 84.66 82.36 90.23 87.05 Mass 95% (%) 80.67 78.95 88.27 85.12 94.46 89.40 Mass Worst (%) 81.23 79.55 90.92 86 31 96.66 92.67 Once again the results between the two filters are similar. Even with highmea s urem e nt noise levels, both filters perfor1ned adequately, missing only one or two times in the Monte Carlo run. The fuel usage numbers are also comparable. Thes e re s ults verify the re s ult s of the pdf compari s on between the filt e r s, as the probability distribution s are similar between th e EKF and SDREF. Therefore, it can be assumed that the filt e r s are
PAGE 97
86 Monte Carlo run. The fuel usage number s are also comparable. These result s verify the results of the pdf comparison between the filters, as the probability distribution s are s imilar between th e EKF and SDREF. Therefore it can be assumed that the filter s are using approximat ely the s ame set of di s tribution functions even though they are achi eve d in different ways. 6.4 4 Initialization Errors Using the new initialization method described in section 6 3.1, a closedloop analysis of the EKF and SDREF was perforn1ed. The initial conditions used for the open loop initialization s tudy were used for thi s inve s tigation The filters were run closedloop in a 3 1 run Monte Carlo simulation. Table 13 display s the re s ults of these runs. Using this new initialization method not only shows better performance in the singlerun, open loop case, but perfo rms equally as well in the closedloop analysis. Probability of hit miss di s tance and f uelu s age number s s how that the SDREF performance is si milar to that of the EKF. Table 13 Miss distance perforrnance for varying initialization error Parameter s (units ) 1 % Error EKF SDREF 5 % Error EKF SDREF PHIT 1 0 1 0 1.0 1.0 CEP ( m ) 0.436 0.377 0 .33 6 0.389 Mean Mis s ( m ) 0 .464 0 .38 0 0.423 0 43 2 95% Miss (m) 0 849 0 .7 50 0.858 0 .77 0 Worst Miss ( m ) 1.01 0 793 1.13 1 .12 Mass CEP (%) 75.13 79.18 72.53 92. 03 Mass Mean (%) 75.00 79.45 72.53 92.03 Mass 95 % (%) 77.34 81.92 74.56 93 .2 0 Mass Wors t (%) 78.12 84.78 75.55
PAGE 98
87 6.5 Summary of Exoatmospheric Guidance Problem The exoatmospheric guidance problem analysis demonstrated, with some minor caveats the SDREF was capable of satisfying the interceptor guidance problem. The assumptions of detectability in solving the SDRE required the addition of a range measurement. An ad hoc initialization algorithm needed to be developed to properly initialize the SDREF With these two restrictions in place however the SDREF performed similarly to the EKF for the exoatmospheric guidance problem for several engagements, noise sources and initialerror combinations.
PAGE 99
CHAPTER 7 PENDULUM PROBLEM FILTER DEVELOPMENT To evaluate the effect of nonlinear dynamics on the SDREF performance, a s imple twodimen s ional pendulum problem was constructed. The scenar io consists of a simple rigid pendulum with friction about the connection point, see (Figure 36). The measurements are assumed to come from an accelerometer attached to the bob, measuring e. 0 0 Figure 36. Pendulum set up 88
PAGE 100
89 The pendulum is given some initial position and velocity and the friction term damps the oscillations. This system, while simple, demonstrates nonlinear dynamics and nonlinear measurements together in the same problem The system dynamics for the pendulum problem are given by the 2nd order nonlinear differential equation . g 0 =sin0B8 L where B0 is the friction ter1n. B was set to be 0.5. Both the SDREF and EKF were derived using the stochastic nonlinear system x=f(x)+rw z = h(x)+v where wand v are Gaussian, zeromean whiteprocess and measurement noise respectively, given by If E[ w(t)wr ( t + r) = W(t)b'(r)] E[ v(t)v(t + r) = R(t)b'( i)] then the state equations for the pendulum system are The measurement equation is z = _.f sin(.x, )Bx2 L (7 .1) (7.2) (7.3) (7.4) (7. 5 ) (7.6)
PAGE 101
90 Therefore, 0 X2 f(x) = g Bx smx L I 2 (7.7) h(x) = g Bx s1nx L I 2 7.1 SDREF Derivation Using the statedependent coefficient parameterization technique described in the exoatmospheric guidance problem, the state equations can be brought to the state dependent linearlike form with 0 K(x) = P(x)H1 (x)R1 F(x) = g sin(x1 ) L X1 and H(x)= B and the resulting filter equations become where i = F(x)x + K(x)[ z(x) H(x)x] K(x) = P(x)Hr (x)R1 and P(i) is the solution of the statedependent Ricatti equation F(x)P + ppT (x)PHT (x)R 1 H(x)P + rrwr = 0 1 B (7.8) (7.9) (7.10) (7 .11) (7 .12) Although this system only has two states and seems simple, a closed for111 so lution of the Ricatti equation prove s difficult. After hours of computation with the symbolic solver Macsyma, a closedloop solution was determined. However the three equation s for .Pi.1 Pi.2 and P2 2 involved over 300 lines of computations. Although these
PAGE 102
91 calculations provided the correct result, it was deterrnined that the previous method used for computing the solution to the guidance problem's Ricatti equation would be used It consisted of a continuous SDRE solver obtained from the LinPak algorithm set available on the worldwideweb. It was easier to incorporate into the pendulum simulation, and executed faster. 7 .2 Extended Kalman Filter Derivation The EKF uses the same stochastic system as the SDREF, equations (7.2) through (7.6). However, the dynamic and measurement matrices are linearized instead of being parameterized as in the SDREF. The linearized equations are given by: and H = oh(x) = ox So the filter equations are 0 g cosx L i g cosx L I I B B i = f (x) + PHT R 1 [ zh(x)] where (7.13) (7 .14 ) (7 .15 ) (7. 16)
PAGE 103
CHAPTERS PENDULUM PROBLEM SIMULATION RESULTS The pendulum problem was analyzed in the visual simulation tool, VisSim The model u sed cons i s t e d of the following parameter s. L = 1.0 (ft) g = 32.2 (ft/ sec2 ) R = 2 (ft/ sec2)i T .05 0 2 r wr = ( rad / s ec ) 0 .05 ( 8 1 ) The simulation was run u s ing fourthorder Runge Kutta integration at a rate of 2 000 hertz ( Hz). The fir s t analy s i s was perfor1ned with both filters receiving the true initial conditi o n s a s their initial estimates The initial condition s used were x1 = 1.0 rad x1 = 1.0 rad X i = 5 0 rad / sec X i = 5.0 rad / sec Figures 37 and 38 show the re s ults of this run. As can be seen in the figure s, both filter s perfor1ned w e ll for thi s scenario Figures 39 and 40 show that the error s for both filter s are almost identical The goal of running this pendulum problem, however was to determine if the pdf function s for the two filters were similar when nonlinear dynamic s were incorporated with the nonlinear meas urements. The PDF curve for each filter, Figures 41 through 46 was determined at three time location s 0.5, 2 0, and 3 0 s econds into the simulation These distributions were computed in the same manner as explained 92
PAGE 104
1.5 .. 1.0 .. . 0.5 'U Ct1 I,... Ct1 0.0 Q.) ..c .. I0.5 .. 1.0 ... 1.5 6 4 . ,. (.) Q.) 2 en 'U .. Ct1 I,... 0 0 .. 0 2 Ct1 Q.) ..c 4 I... 6 ... 8 93 x1 x1 ekf x1sdref I' f\ /", I \ I \ \ I I I 0 1 2 3 4 5 Time (sec) Figure 37 EKF and SD REF estimates vs. true theta with no initialization error I x2 n x2ekf x2sdr I\ A I \ \ I V v I 0 1 2 3 4 5 Time (sec) Figure 38. EKF and SDREF estimates vs. true thetadot with no initialization error. .... f
PAGE 105
0 0 1 0 .. 0.005 ,. 0 0.000 C\1 ,._ C\1 Q) .c. 0 0 0 5 ,.. 0 .010 .. 0.0 1 5 0 .04 0.03 0 .02 u 0.01 Q) en 1J C\1 ,._ 0.00 0 0 0 0 1 C\1 Q) 0.02 .c. 0.03 0.04 ,. .. ,. .. ,. ,. ,. 94 x 1ekf e r r x 1 sd r eerr I I ,, o I ::! ' I I :fl I I I 1iA I I j . l ~ I r ~ < ' I '" u I 1 ; ,, I I w 1tff 0 1 2 3 4 5 Time (sec) Figure 39 EKF and SD REF e s timation error for th e t a with no initiali z ation error x2ekf err x2sdree rr .. I \ I I 'tJ .' I I I / ) I v.: I I y . ' I I \ l I I I I I I 1 I \ \ I '\1', V \. I I ~ ' I \ ,. I I I I I \ ' ., 0 1 2 3 4 5 Time (sec) Figure 40. EKF and SDREF estimation error for theta dot with no initialization error
PAGE 106
8 0 00 6000 ,..._ "l 4 000 V'I 0 e 0. 2000 0 20 1 5 ,..._ V'I 10 V'I 0 e 0. 5 0 .. .. .. .. .. .. 95 EKF ............... SDREF . . . I . 1 . . / / . / 1 / . . .. ... 0 16 0.15 0 14 0 1 3 e (deg) Figure 41 Monte Carlo probability den s ity function for theta with no initialization error time = 0.5 sec. EKF ... ..... ... SDREF . . . . . . . I \ . ,1 .. .. \ / . . . . . t,....._ . .. . 6 .70 6 .68 6 .66 6.6 4 6.62 6 6 0 6 .5 8 6 5 6 6 5 4 e ( deg/ sec ) Figure 42. Mont e C a rl o pro bability den s ity functi o n for theta dot with no initiali z ation e rror time= 0 5 s e c.
PAGE 107
96 3500 3000 EKF ... / '\ SDREF ................ r2500 2000 ... . . . I ...... r . J 1 0 1500 .....\ 1000 ... 500 ... 0 0 .86 10 0 0 I 11/ / . .. . .. . .. 0.8 5 0 .84 0 .83 e ( de g r ees) Figure 43 Monte Carlo probability den s ity function for theta with no initialization error time= 2. 0 sec EKF 0 8 2 / ~ ............... SDREF 800 / .. .. . 600 /: r. 400 I \ 200 I / \ .. V .. ./ 0 I 0.9 6 0.97 0 .98 0 .991. 0 0 1 .0 1 1 .0 2 1 .031.04 1 0 5 1 0 6 1 .0 7 1 .08 1 09 1 1 0 1.11 e ( d eg/ sec) Figure 4 4. M onte Carlo probabilit y d e n si ty func ti o n fo r t heta d o t with no initiali z ation erro r time = 2. 0 sec.
PAGE 108
M N M 0 0 '' 0. N "' d e a. 3500 .. 3000 .. 2500 .. 2000 .. 1500 .. 1000 .. 500 .. 0 0.86 . 1/ ./ . 0.85 ; 97 / '\ . ... . . . . 0 .84 e (deg) . ' EKF ............. SDREF . . .... 0.83 0 .82 Figure 45 Monte Carlo probability den s ity function for theta with no initialization error, time= 3.0 sec. 1000 .. EKF 800 SDREF / ............... ~ . I/. . . .. 600 I: .. 400 I \ 200 .. / I/ 0 .. V .. ..,,...,, I"..._ ' 0 .96 0 .97 0 .98 0 .99 1 .00 1.01 1 .02 1 .03 1 .04 1 .05 1 .0 6 1 .07 1.08 1 .09 1 10 1.11 0 (deg/ sec) Figure 46. M on t e Carlo pro bability density f unction for thetadot, with no initialization error, time= 3.0 sec.
PAGE 109
98 earlier in the guidance problem, Section 6.2. Figures 41 through 46 show, the curves again are similar. The second analysis examined any differences resulting when an initial error was given to each filter The initial conditions used were x1 = 0.0 rad x1 = 0.9 rad x2 = 5.0 rad / sec .x2 = 5.0 rad / sec As was evident in the guidance problem, the SDREF pdf mirrors the EKF pdf once the filter has settled. Figures 47 and 48 show the initial settling differences in the SDREF and EKF. The SDREF is slower to settle, just as in the guidance problem. However, once settled, the performance matches that of the EKF. Once again Figures 49 through 54 demonstrate the similarity in the pdf shape between the two filters. This study help s to confi1n the result s of the guidance analysis. The SDREF, while derived in a different manner leads back to the same internal representation of the pdf as that of the EKF once the initial uncertainty ha s settled.
PAGE 110
0 0 .. ""O Ct1 i.... ..._ 0 5 Ct1 +' Q) .c I.. 1. 0 2.0 (.) 1 5 Q) Ct1 i.... 1 0 +' 0 0 Ct1 +' Q) 0 5 .c I0 0 .. .. .. 99 x1ekferr , ___ .. __ x1sdreerr \ ,'/ ~, ..,., /, I I I t I t I ' I I I I I I I I I I t 0 1 2 3 4 5 Time (sec) Figure 4 7. EKF and SD REF theta estimation error with initialization error x2ekferr x2sdreerr I I I I I I I I I I I I I I ' I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I !(\ ' I I I I I I I I I I I \ \ ' 0 1 2 3 4 5 Time (sec) Figure 48 EKF and SDREF thetadot estimation error with initialization error.
PAGE 111
' 0 N "' 0 e ......... a. ll'l s "1 0 e '" 0. 100 EKF 80 SDREF ......... .. . . 60 . ' . ' . 40 . .. 20 '. . . 0 / I 0 .19 0 .20 0.21 0.22 0.23 0.24 0 .25 0.26 0 .27 0.28 0.29 e ( deg ) 16 14 1 2 10 8 6 4 2 0 .. .. .. .. .. Figure 49. Monte Carlo probability density function for th eta, with initialization error, time= 0 5 sec. EKF ... SOREi . . I I . . . . } \ . 5 05 5.004.954.904 .854.804.754.70 4 .654.604.554.504. 45 e ( deg/ sec) Figure 50 Monte Carlo probability density function for thetadot, with initialization error, time= 0.5 sec
PAGE 112
101 EKF 80 t++++++~ SOREi #. 60 t++i,l,:_. .... _1_ ..,,._ ___ 111~ .. 4 0 t++l++1+++I .. 20t+.... ._.Hl1/'+++_.,,,.\,\... I. +I ,......_ 1/'l s 1/'l 0 e '' 0.. . ~0 ~l1!!~L____::::::::::: I!::".:.:.... _j . 0.560 0 555 0.550 0.545 0.540 0.535 0 .530 0.525 e (deg) 16 14 12 1 0 8 6 4 2 0 Figure 51. Monte Carlo probability den s ity function for theta with initialization e rror time= 2.0 sec. I,_ ,_ .. IEKF SOREi ...... . .Y ...... y \ J \ \ I I I/ \ / 0 .08 0.10 0.12 0.14 0 16 0.18 0.20 0.22 0.24 e ( deg / sec) Figure 52 Monte Carlo probability den s ity function for theta dot, with initialization error, time= 2.0 sec.
PAGE 113
102 100 EKF ... .......... SOREi 80 / ..... '\ ;. I,, ; ' 60 J , ... M ' 0 40 ..._, 0. ...,, "1 V"l 0 0 ..._,, 0. .. 20 \ .. . . ,.: .. ... 0 0.345 0.340 0 335 0 330 0.325 0.320 0.31 5 0.310 16 14 12 10 8 6 4 2 0 ... .. I,, ... .. ... ... ... 0 (deg) Figure 53. Monte Carlo probability density function for theta, initialization error, time= 3 .0 sec. EKF .... . SDREF ,,. '\ :/ \ \ I J \ \ I \ J _,.. V '' 1.52 1.5 0 1.48 1.46 1.44 1. 4 2 1.40 1.38 1.36 1.3 4 0 ( deg/ sec) Figure 54 Monte Carlo probability density function for thetadot, initialization error, time= 3.0 sec.
PAGE 114
CHAPTER9 MSDREFDEVELOPMENT Although the SDREF investigated to this point has demonstrated adequate performance, no proof for stability could be determined. A modification to the fi lter algorithm was found that allowed a stability proof to be formulated. It involved removing the infinitetime horizon assumption in the computation of the covariance matrix. 9.1 MSDREF Derivation Consider the continuous system x = f(x)+Gw = F(x)x+Gw z=h(x)+v = H(x)x+ v (9.1) (9.2) where w is Gaussian zeromean white process noise with E[ w(t)wr (t + i)] = Q(t)8( i) and v is Gaussian zero mean white measurement noise withE[ v(t)vr (t + i)] = R(t)8( i). Also assume that the MSDREF is one where the state update and covariance update are computed simultaneously, as the covariance equation is state dependent Then the equations x = F(x)x + K(x)(zH(x)x) 103 (9.3) (9.4) (9.5)
PAGE 115
104 define the continuou s MSDREF. The system dynamics and filter equations can be rewritten in discrete for1n as with (9.6) ( 9 .7) (9.8) ( 9.9 ) ( 9.10 ) ( 9 .11 ) ( 9 .12) ( 9.13 ) ( 9.14 ) Equations (9.89.12) define the discrete MSDREF. It is also assumed that vk+I and w k+i are random, zero mean inputs and therefore, stability of the MSDREF is independent of these error sources. 9.2 Stability Proof The following proof s hows under the conditions of nonlinear and ''linear' observability, that the MSDREF is locally asymptotically stable. It makes use of a similar proof for the extended Kalman observer (EK0),21 however instead of using the error due to linearization of the EKO, the proof is reformulated to account for H .O.T term s due to a Taylor series expan s ion of the parameterized equation in the MSDREF.
PAGE 116
105 T h eorem: Consider the nonlinear system given in (9.19 2) and equations (9 69.14), if we assume a F k ,Hk are uniformly bounded matrices and Fk i exists ( 9.15) b the nonlinear system is N locally uniformly rank observable; i.e. there exists an integer N 1 such that a rankax hk(x) hk+I (fk (x)) =n ( 9 .16) for all xk EK, where K is a compact set of Rn and hk+t (fk (x)) is defined by the Lie derivative. c the parameterized system is Nlocally uniformly rank observable such that rank Hk (xk) Hk+t Fk (xk) =n for all xk e K, where K is a compact set of R". d. where 1~1fik+1 aik+1 1 + ~1 6k+I for i = 1, ... p and Amax(.) represents the maximum eigenvalue of(.) e. .Jr: ~bjk ~.Jr: for j = 1, ... ,n where ( 9.17) (9.18) (9 .19)
PAGE 117
106 Then the MSDREF i s locally asymptotically stable. Proof: The state estimation error both before and after the update is defined as xk+i = xk+ i xk+t (after update ) xk+ t l k = xk+i xk+llk ( before update ) For ease of notation the following tertns are defined as And finally a candidate Lyapunov function Vk+l is defined V TP, 1k+I = Xk+ l k+I Xk+1 Remark 1 (9 .20 ) (9.2 1 ) (9 22) (9.23) ( 9 .24) (9.2 5 ) For equation (9.25) to be valid, P k must be bounded from above and below A test for the observability of th e nonlinear sys tem22 is a rank dx = n (9.26) x=x.t If n is of full rankfor all x E Rn, then the system is observable for all x. Howev er, satisfying (9.26) d oes not guarantee that the parame terized system used in the MSDREF will be observable. 23 Therefore, t h e observability of the parame t e ri ze d system pair [ F(x),H(x)] must a lso b e test e d Then not only has the nonlinear system been shown to be observable; so that a filter, for which all states are observable, can exist for the
PAGE 118
107 system, but it will also be possible to use the parameterized MSDREF equations. I f both rank co nditions (9.16) and (9.1 7) are met th e co ndition for observability k /Ji / L T ( i, k) H ; T ( RR T ) 1 H ; Cl> ( i k) /3 2 / (9.27) i = k M i s also satisfied, where /31 /32 are positive real numbers, for some finite M 0, andfor all k M, and is bounded from above and below. 24 The condition for which Vk+ i is a decreasing se quence mus t now be deter rnined and will show th e local asymptotic s tability of the MSDREF. The following proof uses2 1 and the form of the MSDREF to establish a theoretical basis for stability. While linearization of the filter equations has proven unsucce ss ful in showing stability of the filter in the past, this approach co nstrain s the state transition and the measurement noi se matrices so as to characterize the error due to parameterization The meas urement residual error, equation (9.24) can be written as (9.28) For each output error prediction co mponent of ek+l' namely e i. k +I' and for each sta te error prediction component of xk+ltk, namely xj,k + l lk; i =I, ... p and j = 1 ... ,n there always exists re s idues r;, k + i and o J,k due to the first order linearization of H i k + i (xk+t) and F1. k (xk) in order to obtain an exact equality for all values of xk+J lk, xk, xk+I' xk, such that e i,k+1 = H i k+1.xk + t 1 k + '1. k + 1 xj,k + 1 t k = FJ,k + 01.k H i ,k+i and Fj,k are the ith and }th rows of H k + t and Fk re s pectively Expanding H i k + i (xk+t )in a Taylor series about xk+i and u s ing the Mean Value Theorem, we have
PAGE 119
108 ei,k+J = H i,k+i (xk+1 )xk+i col' {H1.k+1 (xk+I)} dcolj { Hi.k+I (~ji,k+1)} ,., ,, + a (xk+1 X k+1)} X k+l Xk+l (9.29) j = 1, ... n where the vector ~ji,k+i is that point on the line segment joining the xk+i and xk+i that yields equality in the jth equation of ( 9 29). Let the re s idual error be defined a s then the exact mea s urement error can be defined as or rewritten as where a i,k+J is an unknown timevarying scalar related to r; ,k+ I by 't,k+ i = (1a i.k+i )ei. k+l Through the same process, the equation for the system dynamics i s obtained where j refers to the j 'h component and o j,k is the residual. Written in another equivalent form where bj.k is an unknown time varying scalar related to o j,k by o k = ( b k 1) F kxk } }, }, Recomposing equations (9.32) and (9.35) into vector form yields ak+Jek + i = Hk+1.xk+11k xk+ ltk = b k Fkxk (9.3 0 ) (9.31) ( 9 .32) (9 .33 ) (9.34) (9.35) (9 36) ( 9 .37 )
PAGE 120
109 where ak+l E RPPand bk E Rn.n i s an unknown timevarying diagonal matrix defined as ak+1 = diag(ai, k+1, .. ,ap,k+i ) b k = diag(blk ... ,bnk) If both sides of equation (9.9) are subtracted from xk+i P. H T R I xk+t = xk+Ilk k+t k+1 k+t ek+1 Substituting this into equation (9.27), and using the fact that P k+i = ~ +1 T , ,.., T P. 1T H TR 1 y k+l = xk+l/k k+I xk+llk Xk+llk k+I k+I ek+I T R IH T R IH P. H T R ) ek+I k+I k+lxk+llk + ek+l k+l k+I k+I k+I k+I ek+1 Before proceeding, the inver se relationship i s defined as through a matrix inversion relationship for ~+' ( +) as given in equation (9 .17). Substituting equation (9. 41) into (9. 40), ,, ,_ T{P. 1 H TR 1H }..... T H T R I y k+1 = xk+llk k+1/k + k+I k+I k+I xk+1/k Xk+llk k+1 k+J ek+I T R lH TR IH P. H TR I ek+t k+I k+Ixk+tlk + ek+1 k+1 k+t k+I k+t k+1 ek+t and using the following is obtained V. V ..... TH T R lH ..... ..... T H TR 1 k+l = k+l/k + xk+l/k k+I k+l k+lxk+llk Xk+llk k+I k+I ek+1 T R I H T R IH P. H T R I ek+ I k+l k+lxk+llk + ek+l k+I k+I k+I k+I k+I ek+I Using equation (9.35), equation (9.44) becomes V. V. T ( R I R 1 R 1 k+I = k+llk + ek+I ak+I k+l ak+l ak+I k+l k+I ak+l + R k+t i Hk+t ~+1 Hk+1 r R k+t 1 )ek+i (9.38) (9.39) (9.40) (9.4 1 ) (9.42) (9.43) (9.44) (9.45)
PAGE 121
110 Another way to write Vk+ llk is Using the definition of a decreasing sequence as then V T l T( R I R I R 1 k+l y k = ek+ I ak+I k+I a k + I ak+I k + 1 k+I ak+I + R k+I1 H k + 1 Pk+1 H k + 1 r R k + 1 1 )ek + I +xkr[Fkrbk(FkP k Fkr +Qk)1bkF k P 1k] xk o A sufficient condition to ensure asymptotic stability is that R l R I R I R l H P, H T R l < 0 ak+ I k + l ak+ I ak+ I k + I k + 1 ak+I + k+l k+1 k + I k+1 k+l Conditions which allow this to be true, as s ume that with 1~lfik+I ai,k+l 1 + ~1 fik+l for i = I, ... p (9.4 6 ) (9.47) (9.48) (9. 49 ) (9 50) (9.51) (9.52) where Amax(.) represents the maximum eigenvalue of(.). While the value of ai,k + I cannot be determined a priori the filter can be constructed so that R k + I is chosen such that 0 ~ k + i 1 Similarly It is assumed that F k is invertible and that each bik s atisfies the following condition (9.53)
PAGE 122
111 with (9.54) then equation (9.50) is s atisfied .21 An a priori choice of R k+i and Qk should be chosen to make the interval s of equations (9.51) and (9 .5 3) as large as possible. In the limit ( 9.55) and lim [Jr;, Jr.]= [oo, oo] A, ( QI: )too ( 9 56) The value of a i.k+i can be computed by running the filter and local asymptotic stab ility can then be checked .2 1 Therefore with R k+i and Qk cho se n properly and a i. k+t satisfying equation (9.51) the local asymptotic stability of the MSDREF is shown. No assumption s of a minimum order realization for the estimator are made. It i s obvious that equation (9. 55 ) i s the limiting factor in this approach as the bounds are seve rel y restricted It s hould be noted h oweve r that there is a larg e cla ss of nonlinear physical processes where the meas urement process i s linear therefore yk = Hxk. This make s a i,k+l =IP and equation (9. 55 ) becomes easily sati s fied Q k can then be varied over all possible values to sati sfy equation (9. 56). The refor e a lin ear measurement process allows the local asymptotic s tability conditions to b e met for a much wider set of conditions on R k and Q k
PAGE 123
CHAPTER IO MSDREF PROBLEM SIMULATION RESULTS There are cases in which the SDREF for1nulation has di s tinct advantages over the linearization technique used in the EKF. Consider the sys tem 0 ( 10.1) z = h(x) = X1 +x2 +v (10.2) where x1 ( 0 ) and x2 ( 0 ) are zeromean Gau ss iandi s tributed random variables. Let m and n be the dimensions of the measur e ment and s tate vectors, respectively and h(x) = H(x)x. Before examining the filter performance a proof describing the nonlinear observability of the system is examined. Using a similar methodology,22 and using the complementary estimation terms, this sys tem can be s hown to be weakly observable everywhere except at the origin. Theorem: Given the system equations: and parameterizing as x=f(x) z = h(x) f (x) = F(x)x h(x) = H(x)x 112 (10.3) (10.4)
PAGE 124
113 Let m n. Assume H ( x) has rank n for all x. Then equation ( 10.4 ) is observable for all x, and system given by equation (10.3) is weakly observable on Rn. Proof: The proof is given in21 and the results are used here to demonstrate the observability of the system. If tl1e sufficient condition21 (10.5) is met the system is said to b e weakly observab l e (on S) for all (x ES). U s ing the recursive formula for computi ng Llc ,23 Llc was found to be which can be shown to have rank= l at x1 = x2 = 0. Therefore, the system is weakly observable everywhere except the origin. Three fi lter s will now be used to investigate their performance when confronted with this system. The EKF uses a Taylor series expansion to linearize the measurement matrix about the current state estimate as (10.6) It is easily seen that H(x) lo ses rank whenever x1 = x2 and the system appears unobservable to the EKF algorithm. There i s no sol ution for this problem in the EKF; the SDREF, however has two distinct parameterizations: ( 1 0.7) which can be combined to forrn the parameterized s tate depen dent coefficient
PAGE 125
measurement matrix as H(x,a) = 114 1 A ax2 1 (1a)x1 In this fotm the value of a can be chosen such that los s of rank is avoided. (10.8) This section examines the problem with the MSDREF and EKF. In addition to these two filters a secondorder Gaussian filter ( SOGF) was also derived to extend th e Taylor series expansion of h(x) to include second order ter1ns. As shown in the results, this had very little effect. The equations used for each filter are given below EKF (Continuous): where i = f ( x) + K[ z h( x)] P = F(x)P+ PFr (x) +QPHr (x)R 1H(x)P K= PHr(x)R 1 F(x) = dj(x) dX x=x H(x) = cJh(x) dX x=x SecondOrder Gaussian (Continuous/Discrete): = J(X) +
PAGE 126
where 115 F(x) = aJ(x) ax X X H(x) = ah(x) ax X X Ak = _!_ E[a2(hk ,x()xr ( ) )a2(hk ,x()xr c ) )r J 4 and _!_a2(hk, () )a2(hk, C ) )r 4 ":I 2 a2h ai (h, P ) = tra ce ":I ":11 ox P ax q P. MSDREF (Continuous): i = f (x) + K[zH(x)h(x)] P = F(x)P + ppr (x) + QPHr (x)R _1H(x)P K = PHr (x)R1 ( 10 .12) (10.13) where H(x) and F (x) come from parameterized matrices The three filters were coded and the problem simulated for 10 seconds at a 100Hzupdate rate. Each filter was tuned to give the best perforr11ance that could be achieved. Parameters used for the si mulation runs were When given initial conditions T 0.1 0 V = E[ v( i)v ( r)] = 0 0.01 T 1.0 0 Q = E[w(i)w (r)] = 0 P(O) = 1.0 0 0 1.0 0 0 x(O) = 0 0 x(O ) = 0.1 0.3 1.0 (10. 14 ) (10. 15 )
PAGE 127
116 the EKF and SOGF could not converge to the correct answer. The SDREF, however was able to achieve an accurate estimate of the true states using a = 0.8. This is shown in Figures 55 and 56 The same result s were found for any true initial states where the EKF and SOGF were started with x1 (0) = x2 (0) = 0 which is the logical expected mean value for a zeromean random variable. This caused an algorithmic loss of observability, and the filters converged to the wrong estimate. An example can be constructed where the state estimates never cross the algorithmically unobservable condition of x1 = x2 and therefore the three filters are able to converge properly. One example is 0 2 .x(O) = 0.5 0.1 x(O) = 0.3 Figures 57 and 58 show the adequate convergence for this case. However, the convergence problem should be started at the expected value of the state, x(O) = { O O} r, and achieve convergence when the state estimates cross the algorithmically unobservable condi tion. Only the SD REF f or1r1ulation allows this to be accomplished. This simple academic problem illustrates the usefulness of the SDREF technique when loss of algorithmic observability is an issue. There is nothing that can be done to prevent this in the EKF or SOGF; however the SDREF statedependentcoefficient parameterizations can be adjusted to avoid these conditions.
PAGE 128
X1 117 XHAT O = (0,0) True EKF SOGF SDRE : .. : ~ 02 ,, ,,,. .. ... ;~~~.......,... "':: . ~.,,tJ ~ ,.;,:,:. +.. :. :,;:J .:,:::. '7. :~ ?'7: .. 17'1":::: .. :,:. .. ~,"" . f: .. _,.;. """"9'. ?.; :. ~ ,.:. ;, ,.,;..,~::tt.,i .. '' .. : \!.., ... ,:, ,, ,.r:),,< . ,,, . ',ti .... ,... ... ; .. ,:: ... ', !:.f'' .,, ."' ,..1.. "' .. . ';\f. 'I', :.t: ~; ... t ., tJ "'1,1 r' 1f .. ,t .._ I t 0.0++1++t 0 2 4 6 8 10 Time (sec) Figure 55. True vs. estimated Xl state. 0.4 I I XHAT 0 (0,0) . 0.3 ft.. J \ A A A A A . . . . . .. .. .. .. .. .. .. .. .. ~ vv '\,..~'v.J ''f . f\/"'~vA VV,\ 0.2 .. .. . . .. . . . . ::, .... . ,,.t,.,,,, I \;'". ;.: ' ........ ,:,,;> .. ,. r, .;.::.'\:.; ,,i ':I .. 'r,~ .... ..,. I ..f' \,' ; .,,, ,_ \ ~t"':'l" t :;, ..... t X2 r I 0.1 True . EKF .. SOGF SDRE I 0.0 0.1 0 2 4 6 8 1 0 Time (sec) Figure 56 True vs. estimated X2 state.
PAGE 129
118 0.3 XHAT(0)=(0.2, 0.5) True EKF SOGF SDREF 0.2 .. X2 0.1 ivvl /\,.a . : .. "' ..,, yi1'i.1 ... . . , v V ... .. .. .. . ..,, . ..,,. ,.:_:., ... v X2 . : ' ' 0.0 . 1', 0.1 0 7 2 4 6 8 Time(sec) Figure 57. True vs. estimated X1 state, nonzero initial conditions. 10 0.5... ' ' l XHAT(0)=(0.2, 0.5) l rue EKF SOGF SORE 0 4t+++++~ ' ' < r1~ . 0 .21+.........+.+.~ 0 2 4 6 8 Time (sec) Figt1re 58. True vs. estimated X2 state, nonzero initial conditions. 10
PAGE 130
CHAPTER 11 CRITICAL ANALYSIS OF RESULTS The original purpose of this dissertation was to explore a new nonlinear estimation technique, in hopes of improving upon the limitations of the EKF. Analysis of the new SDRE filtering technique demonstrated an advantage over linearization techniques when confronted with a problem in which loss of algorithmic observability i s encountered. For the other two problems examined, there was little difference in the perf orrnance of the SDREF and EKF filters The first evidence of this occurred with the comparison of the probability distribution functions for the two filters. The SD REF and EKF pdf' s were shown to be almost identical indicating that although the filter models are a1Tived at by, for problems involving full observability, different methods the internal representation of the pdf is similar. This shows that there i s little benefit gained by parameterizing the nonlinearitie s versus using the EKF linearization techniques. This was true even when the filters were exposed to extreme nonlinear conditions for which the EKF might lose track. Again both filters behaved in a similar manner To verify that the use of linear dynamics with nonlinear measurements, did not cause this behavior, the pendulum problem was developed. This allowed examination of the filters when using both nonlinear dynamics and measurements. Again examination of the pdf shape and filter perforrnance yielded similar results. 119
PAGE 131
120 A noticeable shortcoming of the SD REF is the inability to use an initial covariance matrix estimate This proved to be a limitation when the filter is exposed to offnominal initial conditions. It caused an initial oscillation in the state estimates which eventually damped out. This performance, however, was not satisfactory, so an initialization method was created. By assuming that the value of ~, which the SD REF computed is P k ( +),the SDREF can be initialized with an initial covariance () to be used at timestep zero. While this is somewhat of an adhoc technique, it proved successful in removing some of the oscillatory behavior. The oscillatory behavior is due to the fact that the SDREF could be viewed as a multi model adaptive filter, which switches among a bank of steady state filters based on the state estimate .x. When viewed in this way, it is logical to expect that the filter has settled on a steady initial condition When this is not the case, the oscillatory behavior is exhibited. The need for a positive semidefin i te s olution of the Ricatti equation also proved to be limiting factor on the SDREF. While the EKF could handle the bearing s only measurement problem, the SDREF required the addition of a range measurement to satisfy the detectability requirements Once the range measurement was added the filter performance was similar to that of the EKF However, the ability to add measurements to the problem is not always be feasible The SDREF does offer the advantage of multiple parameterizations in hopes of improving performance. This is what allowed the MSDREF to be successful in working the observability problem There is no method to deal with this shortcoming in the EKF. The selection of a good parameterization is similar to the artful ' tuning'' of an EKF. Tuning an EKF requires trial and error tweaking of the process and measurement noi s e
PAGE 132
121 until the desired perfortnance is reached. This would be a similar scenario for the statedependentcoefficient parameterization selection in the SDREF. One of the most important results to emerge from this work is the local asymptotic stability proof for the MSDREF. It was shown that this proof becomes stronger when the system has linear measurements, as the conditions for local stability are easily satisfied for any choice of Rk and Qk. For systems involving nonlinear measurements the bounds on Rk were shown to be somewhat more restrictive. While the new SDREF did not prove to be a substantial improvement over the EKF for most problems, the results of this dissertation proved valuable in understanding why. Through the use of pdf comparisons, sensitivity trades, observability and controllability analyses, and performance studies, the SDREF has been thoroughly investigated. The results are valuable in understanding the internal operations of this filter, its applicability to a real world guidance problem, and the promise of an alternate filtering technique when confronted with loss of algorithmic observability
PAGE 133
REFERENCES 1. Cloutier, J.R. ; Evers, J.H. ; Feeley J.J. ; "Assess ment of AirToAir Missile Guidance and Control Technology," IEEE Control Systems Magazine, Oct., 1989, pp. 2734. 2. Schmidt, G.; "Nonlinear Estimation For Exoatmospheric Trajectories," WL/MN T R 9118 Wright Laboratory, Armament Directorate Eglin AFB FL, Sep., 1991. 3. Gido J.; "Nonlinear Estimation For Exoatmospheric Trajectories : Phase L'' WL/MN TR9109 Wright Laboratory, Armament Directorate, Eglin AFB, FL, May 1991. 4 Gido J.; "Nonlinear Estimation For Exoatmospberic Trajectories: Phase Il, WL/MN TR927003, Wright Laboratory Armament Directorate, Eglin AFB FL, July 1992. 5. Ewing, C.M.; "Bootstrap Estimation: An Application of Sampling Theory to Nonlinear Estimation," Master's thesis, University of Florida, 1992 6. Gordon N.J.; ''Bayesian Method s For Tracking," Dissertation Imperial College, University of London, Feb 1994. 7. Gordon N.J.; Salmond, D.; Ewing, C.M.; Baye s ian State E s timation For Tracking and Guidance Using The Bootstrap Filter," Journal of Guidance, Control, and Dynamics Vol. 18, No. 6, Nov. Dec ., 1995 pp. 14341443 8 Integrated Systems Inc.; '' Intercepting Accelerating Boosters and Reentry Vehicle s Without Range Measurements ," Final Report AFATLTR9099, Oct 1990. 9. Song, T.L .; Speyer, J.L.; ''A Stochastic Analysis of a Modified Gain Extended Kalman Filter with Application to Estimation with Bearing Only Measurements," IEEE Transactions on Automatic Control, Vol AC30, No. 10 (Oct. 1985), pp 940949. 10. Bucy, R.S.; Hecht C.; Senne K D. ; ''An Engineer's Guide To Building Nonlinear Filters ," Vol. 1, SLRTR720004, May 1974 11. Cloutier, J.R.; D 'So uza, C.N.; Mracek C.P.; "Nonlinear Regulator and Nonlinear H00 Control Via The State Dependent Ricatti Equation Technique: Part 1 Theory," Proceedings of the First International Conference on Nonlinear Problems in Aviation and Aerospace, Daytona Beach FL May 1996. 122
PAGE 134
123 12. Cloutier, J.R.; D'Souza, C.N.; Mracek, C.P.; "Nonlinear Regulator and Nonlinear HA Control Via The State Dependent Ricatti Equation Technique: Part 2, Examples," Proceedings of the First International Conference on Nonlinear Problems in Aviation and Aerospace, Daytona Beach, FL, May 1996. 13. Mracek, C.P.; Cloutier, J.R.; D'Souza, C.N.; "A New Technique For Nonlinear Estimation," Proceedings of the 1996 IEEE Conference on Control Applications, Dearborn, Michigan, Sep. 1996 14. Pappano V.; Freidland, B.; ''Speed Sensorless Observer for an Induction Machine ," American Control Conference, June 1997, pp. 38053806. 15. Haessig, D.; Freidland, B.; ''A Method for Simultaneous State and Parameter Estimation in Nonlinear Systems," American Control Conference, June 1997, pp. 947951. 16. Harman, R.; BarItzhack, I.Y.; Azor, R.; ''Satellite Angular Rate Estimation from Vector Measurements,'' Journal of Guidance, Control, a nd Dynamics, Vol 21, No. 3, (May 1998), p. 450. 17. Trailie, J.; "Exoatmospheric Guided Simulation Program," General Electric Company, Philadelphia, PA, March, 1991. 18. Gelb, A.; The Analytic Science Corporation. Applied Optimal Estimation, MIT Press, Cambridge MA, 1986. 19. Stengel, R.F.; Princeton University. Stochastic Optimal Control, John Wiley & Sons New York, 1986. 20. Kwakernaak, H.; Sivan, R.; Linear Optimal Control Systems, John Wiley & Sons New York, 1972. 21. Boutayeb, M.; ''Convergence Analysis of the Extended Kalman Filter as an Observer for NonLinear DiscreteTimeSystems," Conference on Decision and Control, Dec. 1995,pp. 15551560. 22. Isidire, A.; Nonlinear Control Systems, 3rd ed., SpringerVerlag, Berlin, 1995, Chapters 12 23. Hammett, Kelly D.; ''Controllability Issues in Nonlinear StateDependent Ricatti Equation Control," Journal of Guidance, Control, and Dynamics, Vol 21, No. 5 (Se p Oct. 1998) pp. 767773.
PAGE 135
124 24. Deyst, J.; Price C.; ''Condition for Asymptotic Stability of the Discrete Minimum Variance Linear Estimator," IEEE Transactions on Automatic Control Dec 1968 pp. 702705.
PAGE 136
BIOGRAPHICAL SKETCH Craig M. Ewing graduated from The Ohio State University, Co lumbu s, Ohio, in 1986 with a bachelor of science degree in aeronautical/astronautical engineering. He accepted a job with the United States Air Force civil service at Eglin Air Force Base, Florida, where he worked on guidance, navigation, and control issues for the Strategic Defense Initiative for seven years, specializing in the area of guidance and estimation. While working at Eglin, he enrolled at the University of Florida Graduate Engineering Research Center to pursue a master's degree. He received a master of science degree in aerospace engineering in 1992 from the University of Florida. Since 1994 he has been working with Dr. Norman FitzCoy at the University of Florida toward a Ph D. in engineering mechanics. He is presently employed with the simulation and analysis branch of the Munitions Directorate at the U.S. Air Force Research Laboratory, where recent work has focused on the simulation and analysis of new weapon concepts for the Air Force. He is married to the former Darsi Krueger and has two children Tyler, age 5 and Joseph, age 2. 125
PAGE 137
I certify that I have read this study and that in rny opinion it conforms to acceptable sta ndard s of scholarly presentation and is fully adequate, in scope and quality as a disse rtation for the degree of Doctor o '.ft.,..,.J.Lly C Norman G FitzC y, Cha J.J.Pn A ssociate Profe sso r of Aero s ace Engineering Me c hanics, and Engineering S c i e n ce I ce rtify that I have r ea d thi s study and that in my op.in.ion it conforms to acceptable standards of scholarly presentation and i s fully adequate, in sco pe and quality as a disser tation for the degree of Doctor of Louis Cerrato, Jr. Lead Engineer, Aeronautical Systems Center I certify that I have read th.is s tudy and that in my opinion it conforms to acceptable standards of scholarly pre s entation and is fully adequate, in sco pe and quality as a disse rtation for the degree of Doctor of Philo so phy Jame s Cloutier Principal Research Scientist Air Force Re searc h Laboratory I ce rtify that I have read thi s study and that in my opinion i t conforms to acceptabl e s tandard s of scholarly presentation and is fully adequate in scope and quality a s a dissertat ion for the degree of Do c tor of Philosophy ~c. aul A C. Ma s on Ass istant Profes so r of Mechanical Engineering I ce rtify that I have read th.is s tudy and that in my opinion it conforms to acceptable s tand a rd s of scholarly pre se nt a tion and i s fully adequate, in sco pe and qu a lity, as a dissertation for the degree of Doctor of Philo sophy. Pet er H Zipfel Assoc iate Profe sso r of Aero s pace Enginee1ing, M ec hanics, and Engineering Scien ce
PAGE 138
This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. August 1999 C M.J. Ohanian Dean College of Engineering Winfred M. Phillips Dean Graduate School
PAGE 139
3 1262 07056 2292 UNIVERSITY OF FLORIDA II I 1111111 11 11 I I 3 1262 08554 4772

