Estimation of the Orbital States of a Satellite Using the Multivariate Time Series Method

Copyright © IFAC Theory and Application of Digital Control New Delhi, India 1982

ESTIMATION OF THE ORBITAL STATES OF A SATELLITE USING THE MULTIVARIATE TIME SERIES METHOD N. K. Sinha and O. Abul-Haggag Ibrahim Group on Simulation, Optimization and Control Faculty of Engineering, McMaster University, Hamilton, Ontario L8S 4L7, Canada

Abstract: The estimation of the orbital states of a satellite from ground-based meaS'Jrements is compl icated by the fact that a real i stic model of satellite dynamics is nonlinear and so is the observation model. Various methods of nonlinear state estimation have been used in the past. All these methods either require too much computation, or suffer from problems of divergence, In this paper an entirely new approach is presented. This is based on fitting a multivariate time series model directly to the satellite observabIes. This technique has been appl ied to actual data obtained for a satell i te and it has been found that the estimates obtained by this method are very accurate. Keywords: Multivariable systems; modelling; satellites; statistics; orbit determination; forecasting. INTRODUCTION

line estimation purposes. Azim (1979) developed an algorithm combining invariant imbedding with stochastic approximation.

The probl em of estimating the orbi tal states of a satellite is complicated by the fact that a realistic model of satellite dynamics is nonlinear and so is the observation model (Altman, 1972, 1975; Chodas, 1979), This problem has been treated by many nonl inearestimation algorithms (Cox, 1964; Mowery, 1965; Ho and Lee, 1964). Most of these algorithms employ a Taylor series expansion and use the linearized equations to compute the error covariance matrix and the fil ter gains. Sucy (1965) and many others have derived the filtering algorithms, by retaining second-order terms, using the probabilistic approach. Detchmendyand Sridhar (1966) and Kagiwada and others (1969) have derived fil tering algorithms similar to the first-order fil ter for nonlinear estimation problems using the least squares errors cri terion and the invariant imbedd ing technique, These algorithms suffer from a number of major drawbacks, such as the need for prior knowledge of the noise covariance matrices as well as the covariance matrix of the initial estimation error, Also, the major problem usually encountered in nonlinear state estimation is that of divergence.

Recentl y, a swi tching criterion for the combined invariant imbedding-stochastic approx imation algori thm was developed (Sinha and Ibrahim, 1981l. Another approach to the solution of the orbit determination problem is presented in this paper, namely, that of time series analysis. Instead of using the nonlinear dynamic model, a mul tivariate time-series model will be used for forecasting the observables on the basis of the immediate past observations. A stochastic model will be fi t ted to the data obtained from the Communications Research Centre in Ottawa for both slant range and range rate satell ite observ abIes. The model that will be identified and fitted to these two observations is a nonstationary linear stochastic model usuall y referred to as an autoregressive-integrated-moving-average (ARIMA) model, To check the adequacy and lack of fit of the model, several checks will be performed on the residuals to find out if they are in fact whi te noise. Then, one-step ahead forecasts for both series will be obtained, Since the desired orbit is known beforehand, it is possible to use the obtained forecasts to estimate the amount of corrective control action that needs to be taken. Also, in order to obtain state estimates, a Luenberger observer (Sinha and others, 1981) can be used

These problems were the primary motivation for the development of al ternative orbi t determination schemes. The use of secondorder filters was examined (Azim, 1979), but, again these second-order fil ters are not easil y derived for highl y nonl inear systems and most of them cannot be adapted for on-

687

688

N. K. Sinha and O. Abul-Haggag Ibrahim

with the smoothed time-series model.

data

obtained

from

the

rity and invertibility regions, a given spectrum or autocovariance function can be generated by a large number of different ARMA models.

The modelling of the observations by means of a time series model is something less familiar to engineers. Box and Jenkins (1970) introduced a class of linear difference equations called autoregressiveintegrated-moving-average (ARIMA) models for the purpose of univariate stochastic model building. These ARIMA model s can be thought of as representing a time series as the output of a discrete dynamical system to which the input is white noise (MacQ-egor. 1973). Using the B-operator notation. the general mUltivariate ARIMA model of order (p. d. q) may be written as

All these "equivalent" model s (Le .• the ones giving the same autocovariance properties) are said to belong to the same equivalence class. The problem of identification therefore becomes one of devising a set of rules which will select a unique representation from each equivalence class. This means that we have to impose certain constraints on ~ (B) and a (B) and their

MULTIVARIATE TIME SERIES MODEL

~ p (B)

= a q (B)

vd y -t

at

( 1)

-

or. alternatively. (I - ~lB - .. • -~pBP) VdIt

=

(I - a 1B - •.. - aqBq)~t

(2)

where B is the backward shift operator de fined by BIt = I t _l and BSIt = It-s' Also. It E Rn is the observation vector containing the satellite observables at time n

t and {~t ER; t = 1. 2. ...} is a multivariate white noise sequence (or a sequence of random shocks) with mean zero and covariance matrix Ea' ~1' .... ~p' al' .... and aq are matrices of dimensions (nxn) containing the autoregressive and moving average model par ameter s . The general autoregressive term ~ (B) p

d

n

d.

= diag(l-B)

E d. roots i=l 1 on the unit circle to allow the process to be nonstationary (Le .• to allow its level to drift with time). and a stationary autoregressive term ~ (B). Thus. d . is the 1

p

q

different orders and they used sequential identification techniques to estimate these orders. Since the number of parameters can become quite large in these n-dimensional models and because the a-parameters enter in a nonlinear fashion. the parameter estimation step can be very d ifficul t and time consuming. therefore from an estimation viewpoint there is an incentive to keep the order of the moving average terms to a minimum. For this reason. we will restrict ourselves to using pure autoregressive models of the form P

d

(I - tlB - ••. - ~pB )V It

= ~t

(3)

Even though this form may be much less parSimonious than using mixed ARMA canonical forms. the ease of parameter estimation is often argued to override this objection.

vd

consists of two terms; a backward difference term V

p

matrix coefficients which are such that within each equivalence class only one model of this form will satisfy these constraints. Of course. some model structures. or canonical forms. will be more parsimonious than others. Box and Tsao have used full ~p(B) and aq(B) matrices but with elements of

containing

1

degree of differencing necessary to make each observable stationary. The diagonal elements of the matrices ~1' ~2' .... and ~p relate each of the satellite observables to its current and past values. whereas the off-diagonal elements relate each observation to the current and past v al ues of the other observations. To a univariate time series. for a given spectrum or autocovariance function there corresponds a unique ARIMA model providing we restricted ourselves to the stationarity and invertibility region (Box and Jenkins. 1970). However, for a multivariate time series Ciranger and Newbold (1977) showed that even if we restricted ourselves to the stationa-

Identification of Model Orders (p.d) Observations for both slant range and range rate were obtained from real tracking data supplied by the Communications Research Centre in Ottawa. Canada. The sampling rate of the observations is 10 seconds. A sample consisting of 300 observations will be se 1 ec ted and used for build ing the mUltivariate time series model. Now that we have chosen the canonical form which we shall use. the next step is to determine the model orders. p and d. The degree of differencing d. necessary to achieve stationarity is reached when the d.

wi t = V 1 Yit' (i = 1. 2 •...• n) dies out fairly quickly. In a recent paper (Sinha and Ibrahim. 1980). a simple procedure was proposed for determining the degree of differencing most sui table for a given time series. The algorithm used was referred to as the "running average algorithm". where the running average series is obtained from the original series by the following equation autocorrelation function

0

f

Estimation of the Orbital States of a Satellite TABLE 1

k

A. = lk

~

E Yi ' i = 1, 2, •.• , n r= 1 r

689

Choice of Order (p)

(4 )

Model Order (p)

AIC

1

-12018 -12112 -12283 -12313 -12331 -12342 -12341 -12351 -12355

where

A. is the kth running average value lk corresponding to the ith observable Yi. A recursive equation for updating the running average is easily obtained, and is given by

= k+T (kA .

lk

+ y.

)

2

3 4

5 6

(5 )

1

lk+1

8

9 The degree of differencing is then determined from plots 0 f A. vs. k. If the resul ting lk relationship is a curve with zero slope, a constant slope or a constant change in slope, then the degree of differeflcing is zero, one or two, respectively. This algorithm was appiied to both the slant range and range rate series and the degrees of differencing d 1 and d were found to be equal to one and 2 two, respectively. To determine the order of the autoregressive term, p, the fitting process is done by successively fitting autoregressive models of increasing order and evaluating a certain information criterion for each model. This criterion was developed by Akaike in 1916 and is based on information theor y for compar ing different models that have been fitted by the maximum likelihood (ML) method. Our objective is to choose the order p which minimizes the AIC, where AIC = -2 In (maximun of likelihood function) +2 (nunber of estimated parameters)

Estimation of Parameters The next step is to, estimate the values of model parameters, ~ s, as well as of the covariance matrix E of the residual vector ~t. In a manner asimilar to that of the univariate AR-model case, an expression for the autocovariance function of the mul tivariate AR(P)-model can be obtained. This is achieved by multiplying both sides of Eq. (8) by I

'T

_ t j obtain

rj

and

T

= r j_1

~1

In our case the AIC will have the form

(8)

where y is a vector containing the first and second-aifferences of the slant range and the range rate series, respectively.

the

+ ••• + r j _ p

= 1,

r .' s J

we would equations

get

••• , p

2,

by

T

~p'

their

( 1 0)

It It+j

the

(9)

sample

'T

= N t=E1

familiar

Yule-Walker

AT + C.

J-P

= 1,

j

~p'

2, ••• , P

(11)

which together with the equation for j = 0

AT

AT

E = Co - C ~1 - • •• -C ~ a -1 -p P can be solved

,

~2

1 N-j

J

Thus, the order of the autoregressive process which mlnlml zes the AIC is p = 8. The above analysis yields a multivariate model of the form

T

+ r j _2

If we replace estimates

C.

where E is the estimate of the covariance matr ix 8f the resid ual vector ~t' N is the number of observations used and K is the number of estimated parameters. The first term of the criterion aims to selecting the model with minimun residual variance, and the second term penalizes the number of parameters used in fitting the multivariate model. The autoregressive order p and the corresponding values of the AIC are given in Table 1.

taking expectations to

j

(6 )

(1 )

then

to

estimates for the

provid,.e approximate

~'s

the estimated

ML

and Ea.

The results· of parameter summarized in Table 2. Also

(12)

estimation

are

variances of the two A2 residual sequences were est imated as 0a 1 = -6 A2 -13 1.56x10 and 0a2 = 3.88x10 •

690

N. K. Sinha and

TABLE 2 11

j

~j

-.1508 .2188 .3892 .1899 .1505 .1284 .0590 -.0247

2

3 4

5 6 7 8

Parameter Estimates 12

~j

o.

-380.6961 -424.3392 -169.9477 14.7947 21. 0163 -266 . 7457 -256.7389 -282.2958

Check #2

21

22

~j

o. o. o. o. o. O. O. O.

Abul-Haggag Ibrahim

~j

-1. 2682 ·-1.2524 -1. 2269 -1.0588 .8696 - . 6361 - .4448 - .1920

DIAGNOSTIC CHECKS The checks on model adequacy are all aimed at the investigation of the residuals to verify the assumption that these resid ual s have a multivariate GJN distribution with zero mean. Check 111

This check is often referred to as the "Portmanteau lack of fittest", where a Q-criterion is calculated Q

=N

(15 )

If the fitted model

is adequate, then Q is 2

approximately distributed as X (K-p-q). For the first residual sequence, Q was calculated and is equal to 54.81, but from tables 0 f chi-square d istr ibution with K-p-q = 30-8 degrees of freedom and 0.05 level of 2 significance, x.05(22) = 33.9. This indicates a slight lack of fit which again is due to the missing moving average terms. The second residual sequence gave a considerably lower Q2 of 25.79. Thus, there is no evidence of lack of fit. Check 113

The correlation coefficients of the sequence {a} .)} can be evaluated as follows

re i)

=

i N-i

( 13)

E a . (k) a . (k+i) k=1 J J

where r( i) is an estimate of the theoretical covariance of the sequence {a.(.)}, then the estimated autocorrelation is glven by ( 14) These estimated autocorrelations should be uncorrelated and distributed approximately normally about

l

zero

with variance n- , and -112 hence, with a standard error of n These facts could be used to assess approximately the statistical significance of apparent departures of these autocorrelations from zero. Investigation of the auto- and partial correlations of the first residual sequence at showed that the first two autocorrelation 1

values lie beyond the 95~ confidence limits, whereas the partial correlations tend to tail off and damp out exponentially. The obvious reason for this is that we restricted ourselves to a pure autoregressive model form; and the appearance of moving average terms (as indicated by the correlation patterns) was expected. As for the second residual sequence a t ' their corresponding 2

correlation patterns showed no correlation values beyond the 95~ confidence limits. Thus, there seems to be no evidence of model inadequacy from the behaviour of individual correlations.

This final check is concerned with the crosscorrelations between the two residual sequences at and at. This is to validate 1

2

the assumption that the covariance matrix La of the residual sequence has zero-off diagonal elements. The estimated crosscorrelations were all inside the 95~ confidence limits. Also an S-criterion similar to the Q-criterion of Eq. (15) was evaluated using those estimated crosscorrelations yielding S = 24.38, whereas 2 x.05(31-8) = 35.2. From the previous checks we conclude that the mul tivariate model of satellite observables, given by Eq. (8), is fairly adequate. FORECASTING After passing the different diagnostic checks the model is now ready for use in for ecasting. The forecasts that will be obtained are of the minimllll mean square error (MMSE) type. The criterion used for obtaining the forecasts of the observables f-time periods into the future is (16) where It(r) is the estimate of It+f standing at origln t. The above criterion yields the A

MMSE forecast It(r) as the conditional expectation of It+f given knowledge of all the I's up to time t, i.e.

Hence,

for

the mul tivariate

autoregressive

Estimation of the Orbital States of a Satellite

model of Eq. (3), we have ItCf) = t1 I t Cf-1) + t2 I t Cf-2) + ••• + t where

Yt'(f-p)

(18)

p -

1'(.) = ,d 1(.) It(j) = I'(t+j)

= I t ( j)

for j

i

0

otherwise

(19 )

One-step ahead forecasts for the two observables, slant range and range rate, were obtained using Eqs. (18) and (19) with f=1 and t=300, 301, ••• , 324. Figs. 1 and 2 show the actual and predicted observations with the associated prediction errors for slant range and range rate, respectively. CONCLUSIONS The resul ts obtained for forecasting the slant range and the range rate series are very encouraging. It is of importance to point out that the forecast errors resulting from the use of the mul tivariate time series method were significantly lower than those resulting from other nonlinear estimation algoritl'lns. The reason for this is due to the 1 inear i zation process employed in most nonl inear est imation schemes. Another advantage that should not be overlooked is the applicability of the multivariate time series method for the on-line estimation of the orbital observables. This is clear from the fact that the computation time needed per estimate is approximately equal to 0.1 seconds, whereas the sampling rate of the observations is 10 seconds. This gives enough time for the corrective control action to be taken to maintain the satellite in the desired orbit. AI so, the simple structure of the forecast equations allows for on-board microcomputer control, where both the estimation and control steps could be carried out with a small microcomputer installed in the satellite itself. ACKNOWLEDGEMENTS The support of this research by the Department of Communications Canada, under contract OSU19-00056 is gratefully acknowledged. The authors are indebted to Mr. S.P. Altman of the Communications Research Centre for valuable suggestions and to IK. J.F. MacCl'egor of the Department of Chemical Eng ineer ing, McMaster Un i ver si ty, for his helpful suggestions and discussions. REFERENCES Al tman, S. P. (1912). A unified state model of orbital trajectory and attitude dynamics. Ce lestial Mechanics, vol. 6, 425-446. Altman, S.P. (1915). Velocity-space maps and tran:!forms of tracking observations for

691

orbital trajectory state analysis. Celestial Mechanics, vol. 11, 405-428. Azim, S.A. (1919). Nonlinear state estimation with applications to communications satellites. Ph.D . TheSiS, McMaster University, Hamilton, Canada. Box, G.E.P., and G.M. Jenkins (1910). Time Series Analysis-Forecasting and Control: Holden Day, San Francisco. Bucy, R.S. (1965). Nonlinear filtering theory. IEEE Trans. Autom. Control, vol. AC-10, 198-199. Chodas, P. (1919). M.Eng. Thesis, University of Toronto, Toronto, Canada. Cox, H. (1964). CXI the estimation of state variables and parameters for noisy dynamic systems. IEEE Trans. Autom. Control, vol. AC-9, 5-12. Detchmendy, D.M. and R. Sridhar (1966) . Sequential est imation of states and parameters in noisy nonlinear dynamical systems. Trans. ASME, J. Basic Eng . , vol. 88D, 362-366. Granger and Newbold (1911). Forecasting Economic Time Series. Academic Press, New York. Ho, Y.C. and R.C.K. Lee (1964). A Bayesian approach to problems in stochastic estimation and control. IEEE Trans. Autom. Control, vol. AC-9, 333-339. Ibrahim, O.A-H., and N.K. Sinha (1981). A new approach for modelling and forecasting of the slant range for communications satellites. Proc. 12th Annual Pittsburgh Conference on Modelling and Simulation, Pittsburgh, USA, accepted for publication. Kagiwada, H.H., R.E. Kalaba, A. Schumitzky and R. Sridhar (1969). Invariant imbedding and sequential interpolating filters for nonlinear processes. Trans. ASME, J. Basic Eng., vol. 91D, no. 2, 195-200. MacCl'egor, J.F . (1913). Optimal discrete stochastic control theory for process application. The Canadian Journal of Chem. Engrg., vol. 51, 468-418. Mowery, V. o. (1965). Least squares recursive differential correc tion in nonl inear problems. IEEE Trans. Autom. Control, vol. AC-9, 399-401. Sinha, N. K., and O.A-H. Ibrahim (1980). On the choice of the number of integrations in ARIMA models for time series . Prec. Int. Congress on Applied Systems Research and Cybernetics, Mexico, accepted for publ ication . Sinha, N. K. and O.A-H. Ibrahim (1981). The invariant imbedding-stochastic approximation algorithm with application to communications satellites . Proc. 8th Trieannial World IFAC Congress, Japan, accepted for publication. Sinha, N.K . , S. S. Law, and M.H. Li (1981). Microcomputer-based adaptive observer for nonlinear systems. IEEE Trans., vol. IECI-28. Stoica, P. (1911). A test for whiteness. IEEE Trans. Autom. Control, vol. AC-22, No. 6, 992-993.

N. K. Sinha and O. Abul-Haggag Ibrahim

692

39.67792

True Range Estimate

39.67790

:2

39.67788 39.67786 39.67784

o

10

1,f.il'1e igecs.)

Figure l-One-Step Ahear. Forecasts for Slant Panqe

XIO- 6 --- True Range Rate ... Estimate

-2 b 0 • - 280 U Qj

~ -300

:2 - 320 -34 G

-5.

._~

O.OrL-_~ 10

15

__ J 20

Time (secs.) Figure 2-0ne-Step Ahead Forecasts for Range Rate