A method of orbit determination is investigated which employs Picard iteration and Chebyshev series. The method is applied to the problem of determining the orbit of an earth satellite from range and range-rate observations contaminated by noise. It is shown to be readily applicable and to possess linear convergence

Feagin, T.

English

NASA Technical Reports Server

z=.u I
W U)
t$ -j tzi! t 0
* Im CI] n
THE ITERATIVE SOLUTION OF THE PROBLEM -o
OF ORBIT DETERMINATION 4 E o
rb . cH rd 3
tvti
USING CHEBYSHEV SRIES o .
*~OG
-uz H
-rTk -H
FINAL TECHNICAL REPORT Ho e. _
0 (n
tz 0
t4 w *4
Principal Investigator--T. Feagin Ln P OH H
January 1, 1974- January 31, 1975
oc-
Ci)a
University of Tennessee
Space Institute
Tullahoma, Tennessee 37388
NASA Grant NGR-43-OO1-144 -
Sk '-<~
https://ntrs.nasa.gov/search.jsp?R=19750006725 2020-03-23T00:24:02+00:00Z
1. INTRODUCTION
A method of orbit determination is investigated which
employs Picard iteration and Chebyshev series. The method
is applied to the problem of determining the orbit of an
earth satellite from range and range-rate observations
contaminated by noise. The method is shown to be readily
applicable and to possess linear convergence.
2. MATHEMATICAL ANALYSIS OF THE PROBLEM
2.1 FORMULATION OF THE PROBLEM
Consider the mathematical model of a dynamical system
subject to observation given by
x = f(x, t) (1)
= g(x, t) (2)
where
x denotes an n-dimensional state vector;
f, an n-vector function of the state;
y, an m-vector of observations;
g, an m-vector function of state variables; and
t, the time.
The dot denotes differentiation with respect to time.
Let the observations be made at discrete time points
t i (i=0,1,2,. . N) in the interval (t0, tN). Let these
observations be denoted by
(t ), i=0,1,2,. N (3)
In the above model, the functions f and g are nonlinear
in general and the observations are contaminated by measure-
ment errors. The problem of state estimation can now be
stated as follows:
Find the solution x(t) of the differential system (1)
such that
N 2
U = g(x(t ), t i ) - (ti) (4)
i=0
is a minimum. This, in effect, is to find from the possible
solutions of the system (1), the solution that satisfies
the least square criterion with the given observations. If
the system has a unique solution for a given set of initial
conditions, then the problem is to determine the set of
initial conditions, c, for which the solution minimizes the
function U.
2.2 PRESENT METHOD OF SOLUTION
The present work comprises a study of Picard iteration [1]
as a method for the solution of the estimation problem. Picard
iteration has been used to solve the initial value problem [2]
and the two point boundary value problem [3, 4]. In this
case the integration constants are determined such that the
solution satisfies the boundary conditions at each iteration.
The convergence criteria have been established for this pro-
cess. The present work is a natural extension of this method
to the estimation problem. The technique and the method of
implementation are discussed in the following sections.
2.3 PICARD ITERATION
Picard iteration successively approximates the solution
of the differential equation (1) using the following
x i = f(x i-, t) (5)
where i denotes the iteration number and where the initial
(guessed) solution or zeroth approximation is denoted by
x (t). This procedure reduces the solution of the differ-
ential equation to a sequence of simple integrations. It is
also seen that at the end of each iteration the integration
constants c are to be determined. In the case of boundary
value problems, the constants are chosen such that the
given boundary conditions are satisfied by the new approxi-
mate- solution.
For the present problem, the constants are determined
so that the solution xi(t) minimizes the function U given
by Equation 4. With these new values for the constants and
the corresponding new approximate solution, the iteration
process continues to determine the next approximation xi+l(t).
It is expected that the iteration procedure converges
linearly because of its similarity to the classical Picard
iteration.
2.4 MINIMIZATION OF THE RESIDUALS
The minimization of the function U (or the sum of the
squares of the residuals) can be carried out using any of
the conventional schemes to solve the function minimization
problem. For instance, the method of steepest descent could
be used for this purpose [5]. On the other hand, at the
extremum, the function satisfies the following
U N xT 3X q
S = 2 c x [g(xk' t ) - (t ] (6)6c i=O 6c 2x k i
where
6U
- is n x 1 matrix
T
- is an n x n matrix
ac
This is a necessary condition for the extremum of the function
and gives n algebraic equations for the n unknowns in c.
These equations are in general nonlinear in c. However, if
the observations are linear functions of the state, it is
evident from Equation (6) that this equation set is linear
in c (as x is a linear function of c).
The problem of determining c now reduces to the problem
of solving the set of nonlinear algebraic equations (6).
These equations can be solved using Newton's method for
finding the zeros of a function [5, 6]. The derivative matrix
for L = VU can be written as follows:
aLT £ xT gT 6x TgT T T x T
-
-+ (g - Y(t)) (7)
bc 6c xk  )C 6xk  6c bc 6xk ]
This matrix can be used to compute the new value of c using
the following:
-1
j = c - L where (8)
bL
c is an n x n matrix and
L. corresponds to c..J J
2.5 POLYNOMIAL REPRESENTATION
In order to implement Picard iteration, a method of
evaluating the integral (5) must be chosen. The initial
guessed solution, x (t), and the subsequent iterates, xk(t),
can be approximated by a polynomial of k-th degree. A poly-
nomial representation, although other forms are possible,
has the advantage that it can be integrated easily.
For a given degree k, the polynomial Pk(t) of best
approximation to the function c(t) defined on the interval
(-1, 1) minimizes the norm
max [(r) - Pk(T)]  (9)
Polynomials which minimize this norm are said to satisfy the
minimax principle [7]. For many functions, a truncated
series of Chebyshev polynomials is very close to the poly-
nomial of best approximation [7].
For this study, a series of Chebyshev :polynomials. is used to
represent the function xk(t). The Chebyshev polynomial of
k-th degree is given by
Tk() = cos (k arccos (r)) ; -1 r <(+1 (10)
These polynomials can be generated from the recurrence
relations:
Tk+l() = 2 r Tk() k-1l()
To( ) = 1; T (-) T
Picard iteration can easily be implemented using the
orthogonality properties.given in Reference 12. The inter-
val (t , tf) is mapped onto the uniform interval (1, -1) by
the transformation
r = 1 - 2t/(tf - t0 ) (12)
To take advantage of the orthogonality properties, the func-
tions are evaluated at the special points given by
Ti = COS(i w/K) (13)
where K is degree of polynomial representation.
The forcing function f(x i , t) as given in Equation 5
is represented by
K
f(xi, t) = Z b. T. (r) (14)
j=o J
The constants bj can be evaluated by knowing-the function
values of f at the special points T . The integration of
these Equations (14) can easily be done using the integration
relations [7] obtaining the polynomial representation for
xi+1 ()
3. APPLICATION TO ORBIT DETERMINATION PROBLEM
The technique discussed above is employed to solve
the problem of orbit determination with range and range-rate
observations from earth-bound tracking stations. The model
considered and the results are discussed in the following
sections. All computations were performed on an IBM 360/65.
3.1 EQUATIONS OF MOTION
The equations of motion are formulated in an inertial
reference frame. This reference :frame-,,has:: :its origin at,,:the
center of the earth. The X-Y plane coincides with the equa-
torial plane of the earth with the X-axis pointing to the
first point of Aries. The Z-axis is in the direction of the
north pole and the Y-axis forms a right-handed triad with
the other two axes (see Figure 1).
The equations of motion of a satellite, considered
to be a mass point moving in the central gravity field of
the earth, in the inertial system of co-ordinates are given
by [8]
x--3X = -X/R (15)
3
Y= -Y/R (16
Z = -pZ/R (17)
where
R = (X2 + 2 + Z2 ) 1 / 2
4 = gravitational constant (3.986 x 105 km3/sec 2
zNorth Pole
nS
!i P
TS --ith Tracking
Station
S--Satellite
P--Perigee
R--Range
n--Unit Normal to
Orbital Plane
Figure 1. Relative positions of satellite and tracking
station.
3.2 OBSERVATION MODEL
The earth is modeled as a sphere rotating with an
angular velocity of 360 deg/day about the Z-axis. The loca-
tion of the tracking station is specified by giving the
latitude and the longitude of the station measured with
respect to the inertial co-ordinate system (XYZ) as shown
in Figure 1 at the instant t o .
Two typesvof'observations are made by each station:
(1) the range (the distance between the tracking station and
the satellite) and (2) the range-rate (the time rate of the
range at the time of observation) [9, 10]. In a practical
case the range and range-rate are obtained by the radar
network. These observations are denoted by
A^ R i j = 1, 2, . . . S(18)( (18)Y i . (3) i = 0, 1, 2, . . N
Ri
Here, the superscript j refers to the station number with
the total number of stations considered being S.
The functional form of the observation vector g is
yet to be given. Let X (t), Y() (t) and Z ( ) (t) be the
co-ordinates of the j-th tracking station at time t. Then
the observation vector as a function of the state variables
can immediately be written as
R 01(t)
g(t) = (j) (19)
where
R(J)(t) = {X(t) - X( j ) (t)j2 +Y(t) - (j)(t) 2
(20)
+ [Z(t) - Z(J)(t)J2 1/ 2
RCj (t) = ; *ILX(t) -X (t)] [X(t) - X'(t)]
+ [Y(t) - Y()(t)] [Y(t) - Y(j)(t)] (21)
+ [Z(t) - z(J)(t)] [z(t) - z(J(t)]}
From the above it is seen that both the observations
and the equations of motion are nonlinear.
3.3 DETAILS OF NUMERICAL SOLUTION
The Equations 15 to 17 representing the system dynamics
are second-order nonlinear differential equations. These
can be written as a set of six first order nonlinear equations,
one for each of the state variables (viz, X, Y, Z, X, Y, Z)
in the form given by Equation 1. The state variables are
represented by a series of Chebyshev polynomials with 16
terms, the number of terms being chosen based on the recom-
mendations of the earlier studies [1, 11] on the representation
of the orbits using Chebyshev polynomials. The computation
of the forces is performed only at the selected points as
given by Equations 12 and 13 in order to maintain the ortho-
gonality property of the polynomials. The integration of
the second-order system gives rise to two vector constants
of integration (each with three elements). Let these six
constants be represented by c. The function to be minimized
for the problem is given by
U= S N (R (t) - R. )]2 ( (j) ) 2} (22) = {[R (t) -  i I + [R (ti) - R ] 
j=l i=0 i O
The function given by Equation 22 is minimized with respect
to c by finding the solution of the system of equations
given by
6UVU = p = 0 (23)
Newton's method is used to find the zeros of VU. No major
difficulties are encountered in implementing this scheme.
It is observed that U decreases as the solution is
approached.
3.4 NUMERICAL RESULTS AND DISCUSSIONS
To demonstrate the validity of the method and the
scheme of implementation two examples are considered. The
first one is the determination of the state of a satellite
moving in a circular orbit. The orbital elements are given
in Table I. The observations are made in the interval
(0, 1513.1) seconds. This corresponds to one-fourth of a
revolution in the orbit. The tracking station data (viz,
the location and the number of observations) are also pre-
sented in Table I. It should be noted that the observations
are not simultaneous except for the first. This does not in
any way constrain the applicability of this method for a
general set of observations. The observations for both
examples presented here are simulated on the computer using
th
a random number generator. At the i iteration the error is
E = [(X i - X) 2  (Yi - 2 + (Z Z)2 1/2 (24)
The initial guessed solution for this problem is in
error by about 72 km throughout the interval considered. In
Figure 2, the error as a function of time is plotted for the
first two iterations. From this figure it can be seen that
although the initial guess is far from the solution, within
two iterations the error is reduced to less than 5 km. To
have a clear pictorial representation of the convergence of
the process, the logarithm of the absolute error is pre-
sented in Figure 3, as a function of the iteration number.
The error at selected time points is plotted in order to
TABLE I
STATION AND ORBIT DATA (EXAMPLE I)
Station Number
Parameters 1 2 3
Latitude (deg) 18.0 12.0 10.0
Longitude (deg) 0.0 28.0 14.0
Number of observations 10 20 30
Interval between observations
(sec) 168.0 79.0 52.0
Semi-major axis (km) a = 7178.145
Eccentricity e = 0.0
Inclination (deg) i = 20
Longitude of ascending node 2 = 0.0
Argument of perigee a = 0.0
80-
Initial Guess
60-
40-
0
1st Iteration
20-
2nd Iteration
0-
0 500 1000 1500
Time (sec)
Figure 2. Error as. a functioniof time .for different
iterations (circular orbit).
6-
5
4
3-
2
o 1-
4-
w T=1513.1 (sec)
T=1046 1
44461
O
o -2 -
T=0 o 0
-3
-4
T=467.0
-5
1 2 3 4 5 6 7 8. 9
Iteration Number
Figure 3. Convergence of the method for a circular orbit.
depict the convergence of the over-all solution. It is
apparent from this that the error reduction is linear as the
solution approaches the true solution.
The convergence criterion is that the error be less
than a meter. The iteration is terminated when this cri-
terion is satisfied. -For this problem, the convergence of
the process is attained at the end of nine iterations. The
time for the computation is 21.6 seconds.
The second problem considered is the case of a
satellite moving in an elliptic orbit. The orbit has an
eccentricity of 0.0557 with a perigee height of 400 km.
Here, the arc considered for the solution is from t = 0
(perigee) to tf = 1513.1 seconds. The circular orbit of the
first example and this ecdentric orbit have the same semi-
major axes, and therefore, the same orbital periods. This
example allows a study of the effect, if any, of eccentricity
on the method of solution. The tracking stations and the
times of observation are unaltered from those of the first
problem. The tracking station data and the orbit data are
presented in Table II.
Figure 4 contains the graphical representation of the
errors in the solution for the first two iterations. The
error of about 70 km for the guessed solution is reduced
to about 10 km in two iterations. The logarithm of the
absolute error is plotted in Figure 5 as a function of the
iteration number. Here again, it is evident that the con-
vergence is linear. The process converges, for this problem,
in 11 iterations to a solution which is in error by about a
meter.
In summary, it is Eseen that the method is applicable
to the problem of orbit determination and that the conver-
gence is linear as expected. The comparison of the two
problems considered here is presented in Table III. It is
seen that the elliptic orbit takes more iterations.
TABLE II
STATION AND ORBIT DATA (EXAMPLE II)
Station Number
Parameters 1 , 2 3
Latitude (deg) 18.0 12.0 10.0
Longitude (deg) 0.0 28.0 14.0
Number ,of observations 10 20 30
Interval between observations
(sec) 168.0 79.0 52.0
Semi-major axis (km) a = 7178.145
Eccentricity e = 0.0557
Inclination (deg) i = 20
Longitude of ascending node 0 = 0.0
Argument of perigee W = 0.0
40-
20
2nd Iteratio
0 500 1000 1500
Time (sec)
Figure 4. Error as a function of time for different
iterations (elliptic orbit).
5-
4
3-
2-
T=1513.1 (sec)
o 1
ST=1046.1
d4
0 \ T=O
0
-3
-4 T=467 0
-5
1 2 3 4 5 6 7 8 9
Iteration Number
Figure 5. Convergence of the method for an elliptic orbit.
TABLE III
COMPARISON OF THE TWO EXAMPLES
Example I Example II
Parameter . (Circular Orbit) (Elliptic Orbit)
Arc length in time (sec) 1513.1 1513.1
Number of observations 50 50
Eccentricity 0.0 0.0557
Initial error in state
variables 1.0% 1.0%
Number of iterations for
donvergence 9 11
Computer time (sec) 21.6 28.9
4. CONCLUSIONS
An iterative scheme using Picard iteration has been
investigated for the solution of the state estimation pro-
blem with discrete observations. The applicability of this
scheme to practical problems has been demonstrated [12]
using as an example the problem of orbit determination of
an earth satellite with range and range-rate observations
from earth-bound tracking stations.
Unlike some of the more commonly used-methods, this
method does not require the formulation or the solution of
the linear perturbation equations. From the examples con-
sidered, it is seen the method has linear convergence.
This scheme can be extended easily to a general problem
or orbit determination (including the effects of drag, oblate-
ness of the earth, gravity fields of other bodies, etc.).
However, further work is necessary to study the convergence
of Chebyshev series, and the number of terms needed to
represent the solution over a given interval.
A sequential algorithm for estimating the orbit using
Chebyshev series is presently being investigated for use in
a real time orbit determination environment.
LIST OF REFERENCES
LIST OF REFERENCES
1. Falb, .Peter L., and Jaul de John. Some...Successive
Approximation Methods in Control and Oscillation
Theory. New York: Academic Press, 1969.
2. Feagin, T., "The Numerical Solution of Two Point
Boundary Value Problems Using Chebyshev Series."
Ph.D. dissertation, University of Texas, Austin,
1972.
3. Vulikh, B. Z., Introduction to Functional Analysis for
Scientists and Technologists. New York: Pergamon
Press, 1963.
4. Picard, E., "Sur L'application des Methodes
D'approximations Successives a L'etude de Certaines
Equations Differentielles Ordinaires," Journal of
Mathematics, 9:217-271, 1893.
5. Pierre, Donald A., Optimization Theory with Applications.
New York: John Wiley and Sons, Inc., 1969.
6. Hamming, R.. W., Numerical Methods for Scientists and
Engineers. New York: McGraw-Hill Book Company,
Inc., 1962.
7. Fox, L., and I. B. Parker, Chebyshev Polynomials in
:Numerical Analysis. New York: Oxford University
Press, 1968.
8. Baker, R. M. L., Jr., Astrodynamics-Applications and
Advanced Topics. New York: Academic Press, 1967.
9. Deutsch, Ralph, Orbital Dynamics of Space Vehicles.
Englewood Cliffs, New Jersey: Prentice Hall, Inc.:
1963.
10. Escobal, P. R., Methods of Orbit Determination. New
York: John Wiley and Sons, Inc., 1965.
11. Wright, K., "Chebyshev Collacation Methods for Ordinary
Differential Equations," Computer Journal,
6:358-365, January, 1964.
12. Mikkilineni, P. and T. Feagin, "The Determination of
Orbits .Using Picard Iteration," NASA Special
Publication (in print), 1975.


The iterative solution of the problem of orbit determination using Chebyshev series

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server