# Comparison between Two Methods to Calculate the Transition Matrix of Orbit Motion

Authors: Ana Paula Marins Chiaradia; Hélio Koiti Kuga; Antonio Fernando Bertachini de Almeida Prado
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

Two methods to evaluate the state transition matrix are implemented and analyzed to verify the computational cost and the accuracy of both methods. This evaluation represents one of the highest computational costs on the artificial satellite orbit determination task. The first method is an approximation of the Keplerian motion, providing an analytical solution which is then calculated numerically by solving Kepler's equation. The second one is a local numerical approximation that includes the effect ofJ2. The analysis is performed comparing these two methods with a reference generated by a numerical integrator. For small intervals of time (1 to 10 s) and when one needs more accuracy, it is recommended to use the second method, since the CPU time does not excessively overload the computer during the orbit determination procedure. For larger intervals of time and when one expects more stability on the calculation, it is recommended to use the first method.


---

## Body

## 1. Introduction

The orbit determination process consists of obtaining values of the parameters which completely specify the motion of an orbiting body, like artificial satellites, based on a set of observations of the body. It involves nonlinear dynamical and nonlinear measurement systems, which depends on the tracking system, and estimation technique (e.g., the Kalman filtering or least squares [1–6]). The dynamical system model consists of the description for the dynamics of the orbital motion of a satellite, measurement models, Earth’s rotation effects, and perturbation models. According to Montenbruck and Gill [7], these models depend on, besides the state variables that define the initial conditions, a variety of parameters that either affect the dynamical motion or the measurement process. Due to the complexity of the applied models, it is hardly possible to have a direct solution for any of these parameters from a given set of observations. To linearize the relation between the observables and the independent parameters, one should obtain simplified expressions that can be handled more easily. A linearization of the trajectory and measurement model requires a large number of partial derivatives, among them, the state transition matrix [7–9].The function of the transition matrix is to relate the rectangular coordinate variations for the timestk and tk+1. The evaluation of the state transition matrix presents one of the highest computational costs on the artificial satellite orbit determination procedure, because it requires the evaluation of the Jacobian matrix and the integration of the current variational equations. This matrix can pose cumbersome analytical expressions when using a complex force model [10]. Binning [11] suggests one method to avoid the problem of the high computational cost and extended analytical expressions of the transition matrix. This method consists of propagating the state vector using complete force model and, then, to compute the transition matrix using a simplified force model. For sure, this simplification provides some impact on the accuracy of the orbit determination. The analytical calculation of the transition matrix of the Keplerian motion is a reasonable approximation when only short time intervals of the observations and reference instant are involved. On the other hand, the inclusion of the Earth-flattening (J2) effect in the transition matrix can be performed adopting Markley’s method [12].In Chiaradia [13], Chiaradia et al. [14], and Gomes et al. [15], simplified and compact algorithms with low computational cost are developed for artificial satellite orbit determination in real time and onboard, using the global positioning system (GPS). The state vector, composed of the position, velocity, bias, drift, and drift rate of the GPS receiver clock, has been estimated by the extended Kalman filter in all three works. The fourth-order Runge-Kutta numerical integrator has been used to integrate the state vector. The equations of motion have considered only the perturbations due to the geopotential. The state error covariance matrix has been propagated through the transition matrix, which has been calculated considering the pure Keplerian motion. To improve the accuracy of those algorithms, two methods have been compared for calculating the transition matrix considering circular (1000 km of altitude) and elliptical (Molniya) orbits. Those methods considered the pure Keplerian motion and perturbed only by the flattening effect of the Earth. Markley's method is used to include the flattening of the Earth in the transition matrix. It is a method that allows the inclusion of more perturbations in a simple way.In the present work, the spherical harmonic coefficients of degree and order up to 50 and the drag effect are included in the reference orbit provided by the RK78 numerical integrator (Runge-Kutta with Fehlberg coefficients of order 7-8). For this study, the atmospheric density is considered constant. For the simulations made here, circular and elliptical low (up to 300 km) orbits were used. To analyze the results, the reference orbit is compared with the two methods implemented in this work.


## 2. State Transition Matrix

The differential equation for the Keplerian motion is expressed by(2.1)r̈=-μrr3,
with r=(x,y,z),(2.2)r=x2+y2+z2,
where r is the magnitude of the satellite position vector and μ is the Earth gravitational constant. The equation that relates the state deviations in different times is given by(2.3)(δrδv)=Φ(t,t0)(δr0δv0),
where r and v are the position and velocity vectors at the time t, respectively, r0 and v0 are the initial position and velocity vectors at the time t0, respectively, and Φ is the transition matrix given by(2.4)Φ=(Φ11Φ12Φ21Φ22)=(∂r∂r0∂r∂v0∂v∂r0∂v∂v0).
The submatrices Φ11, Φ12, Φ21, and Φ22 are calculated in accordance with the two methods that will be shown in the next topics.


## 3. First Method: The Analytical Transition Matrix Solution Considering the Pure Keplerian Motion

Goodyear [16] published a method for the analytical calculation of a transition matrix for the two-body problem. This method is valid for any kind of orbit. Kuga [10] implemented this method using the same elegant and adequate formulation optimized for the Keplerian elliptical orbit problem. He performed some simplifications in this method to increase its numerical efficiency (processing time, memory, and accuracy). However, this method is used to propagate the position and velocity of the satellite and can be used in any kind of two-body orbits. Such equations are simple and easy to be developed.According to Goodyear [16] and Shepperd [17], the four submatrices 3×3 are obtained developing (2.4). Then, they are written as(3.1)Φ11=fI+(rv)(M21M22M31M32)(r0v0)T,Φ12=gI+(rv)(M22M23M32M33)(r0v0)T,Φ21=ḟI-(rv)(M11M12M21M22)(r0v0)T,Φ22=ġI-(rv)(M12M13M22M23)(r0v0)T,
where I is the identity matrix 3×3, (rv) is matrix 3×2, and the Mi,j, i,j=1,2,3, are the components of a 3×3 matrix M, which will be shown later, and f, g, ḟ, ġ are calculated in the next topic [8, 9].


### 3.1. Calculating the Functionsf, g, ḟ, ġ

Givenr0=(x0,y0,z0), v0=(ẋ0,ẏ0,ż0), and the propagation interval Δt=t-t0,(3.2)r0=x02+y02+z02,h0=x0ẋ0+y0ẏ0+z0ż0,v0=ẋ02+ẏ02+ż02,α=v02-2μr0,1a=-αμ,
where a is the semimajor axis. The eccentric anomaly E0 and the eccentricity e for the initial orbit are calculated by(3.3)esin⁡E0=h0μa,ecos⁡E0=1-r0a,
with E0 reduced to the interval 0 to 2π. The mean anomalies for the initial and propagated orbits, M0 and M, are given by(3.4)M0=E0-esin⁡E0,n=μa3,M=nΔt+M0,
with M0 and M reduced to the interval 0 to 2π. The eccentric anomaly for the propagated orbit is calculated by Kepler’s equation. Kepler’s equation is solved from an initial guess based on an approximation series, then iterated by Newton-Raphson’s method, until convergence to the level of 10−12 is achieved. The variation of the eccentric anomaly is calculated and reduced to the interval 0 to 2π:(3.5)ΔE=E-E0.
The transcendental functions s0, s1, s2 for the elliptical orbit, according to Goodyear [16], are calculated by(3.6)s0=cos⁡ΔE,s1=aμsin⁡ΔE,s2=aμ(1-s0).
Therefore, the functions f, g, ḟ,ġ are valid only for elliptical orbits and are calculated by(3.7)r=r0s0+h0s1+μs2,(3.8)f=1-μs2r0,(3.9)g=r0s1+h0s2,(3.10)ḟ=-μs1rr0,(3.11)ġ=1-μs2r.
The propagated vectors r e v are given by(3.12)r=r0f+v0g,v=r0ḟ+v0ġ.
There is no singularity problem and Kepler’s equation is solved through Newton Raphson’s method in double precision. The classical parameters a, e, i, Ω, ω, are constant in the Keplerian motion, and, therefore, there is no different subscript for them.


### 3.2. The Evaluation of the MatrixM3×3

First of all, it is necessary to calculate the secular componentU including the effect of multi-revolutions in the case where the orbit propagation time, Δt, is larger than one orbital period:(3.13)ΔE=ΔE+INT(Δtn2π)2π,s4′=cos⁡ΔE-1,s5′=sin⁡ΔE-ΔE,U=s2Δt+(aμ)5(ΔEs4′-3s5′),
where INT(x) provides the truncated integer number of the real argument x, and the variables s4′ and s5′ are related with the transcendental functions s4 and s5 (Goodyear [16]). Therefore, the components of the matrix M are related to the partial derivatives of (3.8)–(3.11) with respect to the r0 and v0, given by(3.14)M11=(s0rr0+1r02+1r2)ḟ-μ2Ur3r03,M12=(ḟs1r+(ġ-1)r2),M13=(ġ-1)s1r-μUr3,M21=-(ḟs1r0+(f-1)r02),M22=-ḟs2,M23=-(ġ-1)s2,M31=(f-1)s1r0-μUr03,M32=(f-1)s2,M33=gs2-U.


### 3.3. The Property of the Transition Matrix

Sometimes the inverse matrix is required, such as in backward filters. It is also easily accomplished as follows. The inverse matrixΦ-1, which propagates deviations backward from t to t0, is given by(3.15)Φ-1=(Φ22T-Φ12T-Φ21TΦ11T),
which results from the canonic nature of the original equations [18, 19]. This also applies to the second method (Markley’s) if the perturbations are derived from a potential (e.g., J2 effect).


## 3.1. Calculating the Functionsf, g, ḟ, ġ

Givenr0=(x0,y0,z0), v0=(ẋ0,ẏ0,ż0), and the propagation interval Δt=t-t0,(3.2)r0=x02+y02+z02,h0=x0ẋ0+y0ẏ0+z0ż0,v0=ẋ02+ẏ02+ż02,α=v02-2μr0,1a=-αμ,
where a is the semimajor axis. The eccentric anomaly E0 and the eccentricity e for the initial orbit are calculated by(3.3)esin⁡E0=h0μa,ecos⁡E0=1-r0a,
with E0 reduced to the interval 0 to 2π. The mean anomalies for the initial and propagated orbits, M0 and M, are given by(3.4)M0=E0-esin⁡E0,n=μa3,M=nΔt+M0,
with M0 and M reduced to the interval 0 to 2π. The eccentric anomaly for the propagated orbit is calculated by Kepler’s equation. Kepler’s equation is solved from an initial guess based on an approximation series, then iterated by Newton-Raphson’s method, until convergence to the level of 10−12 is achieved. The variation of the eccentric anomaly is calculated and reduced to the interval 0 to 2π:(3.5)ΔE=E-E0.
The transcendental functions s0, s1, s2 for the elliptical orbit, according to Goodyear [16], are calculated by(3.6)s0=cos⁡ΔE,s1=aμsin⁡ΔE,s2=aμ(1-s0).
Therefore, the functions f, g, ḟ,ġ are valid only for elliptical orbits and are calculated by(3.7)r=r0s0+h0s1+μs2,(3.8)f=1-μs2r0,(3.9)g=r0s1+h0s2,(3.10)ḟ=-μs1rr0,(3.11)ġ=1-μs2r.
The propagated vectors r e v are given by(3.12)r=r0f+v0g,v=r0ḟ+v0ġ.
There is no singularity problem and Kepler’s equation is solved through Newton Raphson’s method in double precision. The classical parameters a, e, i, Ω, ω, are constant in the Keplerian motion, and, therefore, there is no different subscript for them.


## 3.2. The Evaluation of the MatrixM3×3

First of all, it is necessary to calculate the secular componentU including the effect of multi-revolutions in the case where the orbit propagation time, Δt, is larger than one orbital period:(3.13)ΔE=ΔE+INT(Δtn2π)2π,s4′=cos⁡ΔE-1,s5′=sin⁡ΔE-ΔE,U=s2Δt+(aμ)5(ΔEs4′-3s5′),
where INT(x) provides the truncated integer number of the real argument x, and the variables s4′ and s5′ are related with the transcendental functions s4 and s5 (Goodyear [16]). Therefore, the components of the matrix M are related to the partial derivatives of (3.8)–(3.11) with respect to the r0 and v0, given by(3.14)M11=(s0rr0+1r02+1r2)ḟ-μ2Ur3r03,M12=(ḟs1r+(ġ-1)r2),M13=(ġ-1)s1r-μUr3,M21=-(ḟs1r0+(f-1)r02),M22=-ḟs2,M23=-(ġ-1)s2,M31=(f-1)s1r0-μUr03,M32=(f-1)s2,M33=gs2-U.


## 3.3. The Property of the Transition Matrix

Sometimes the inverse matrix is required, such as in backward filters. It is also easily accomplished as follows. The inverse matrixΦ-1, which propagates deviations backward from t to t0, is given by(3.15)Φ-1=(Φ22T-Φ12T-Φ21TΦ11T),
which results from the canonic nature of the original equations [18, 19]. This also applies to the second method (Markley’s) if the perturbations are derived from a potential (e.g., J2 effect).


## 4. Second Method: Markley’s Method

Markley's method uses two states, one at thetk-1 time and the other at the tk time, and calculates the transition matrix between them by using μ, J2, Δt, the radius of the Earth, and the two states. In this case, the effect of the Earth’s flattening is the most influent factor in the process. Markley’s method consists of making one approximation for the transition matrix of the state vector based on the Taylor series expansion for short intervals of propagation, Δt. This method can be used by any kind of orbit and the equations are simple and easily implemented, as shown next. The state transition's differential equation is defined by (4.1)dΦ(t,t0)dt=A1(t)Φ(t,t0)=[0IG(t)0]Φ(t,t0),
where Φ(t0,t0)≡I is the initial condition, r=(xyz)T and v=(ẋẏż)T are the Cartesian state at the instant t, r0 and v0 are the Cartesian state at the instant t0, 0≡ the matrix 3×3 of zeros, I≡ the identity matrix 3×3,  G(t)≡∂a(r,t)/∂r≡the gradient matrix, and a(r,t)= the accelerations of the satellite.Performing successive derivatives of (4.1), followed by substitutions, gives the derivative of the transition matrix:(4.2)diΦdti=Ai(t)Φ(t,t0),
where
(4.3)Ai(t)=Ȧi-1(t)+Ai-1(t)A1(t).
The dot represents the derivative with respect to the time. Developing Φ(t,t0) in Taylor's series at t=t0, using the matrices Ai(t0) for i=1,…,4 and the initial condition Φ(t0,t0)≡I, the transition matrix of the position and velocity obtained after some simplifications is given by [12](4.4)Φ(t,t0)≈[ΦrrΦrvΦvrΦvv]6×6,
where(4.5)Φrr≡I+(2G0+G)(Δt)26,Φrv≡IΔt+(G0+G)(Δt)312,Φvr≡(G0+G)(Δt)2,Φvv≡I+(G0+2G)(Δt)26,Δt≡t-t0,G0≡G(t0).
The calculations of these matrices pose no problems, since the gradient matrix G in the end of the propagating interval is a function of the final Cartesian state which should be calculated during the data processing and is available at no extra cost. The G and, therefore, Φrr,Φrv,Φvr,Φvv are symmetric if the perturbation derives from a potential. The G gradient matrix including only the central force and the J2 is given by(4.6)G(t)=∂a(r,t)∂r=[∂ax∂x∂ax∂y∂ax∂z∂ay∂x∂ay∂y∂ay∂z∂az∂x∂az∂y∂az∂z].
The accelerations due the central force and the Earth’s flattening are given by(4.7)ax=-μxr3[1+32J2re2r2(1-5z2r2)],ay=yxax,az=-μzr3[1+32J2re2r2(3-5z2r2)].
The partial derivatives are [20](4.8)∂ax∂x=μr5[3x2-r2-32J2re2+152J2re2r2(x2+z2)-1052J2re2r4x2z2],∂ax∂y=3μxyr5[1+52J2re2r2-352J2re2r4z2],∂ax∂z=3μxzr5[1+152J2re2r2-352J2re2r4z2],∂ay∂x=∂ax∂y,∂ay∂y=yx∂ax∂y+axx,∂ay∂z=yx∂ax∂z,∂az∂x=∂ax∂z,∂az∂y=∂ay∂z,∂az∂z=μr5[-r2+3(z2-32J2re2+15J2re2r2z2-352J2re2r4z4)].


## 5. Analyses of the Methods

First of all, the reference orbit is integrated using the numerical integrator Runge-Kutta of eighth order with automatic step size control (RK78). This integrator is implemented to integrate simultaneously the stated vectors considering the spherical harmonic coefficients up to 50th order and degree and the transition matrix considering the Earth's flattening and the atmospheric drag effects.The transition matrix generated by the numerical integrator is used as the reference to compare the transition matrix generated by the two methods. To compare the accuracy of them, let us define the global relative error as [10](5.1)εGlobal=136∑i=16∑j=16|(Φij-Φ̅ij)Φ̅ij|,
where Φ̅ij is the component i, j of the transition matrix calculated by the RK78 considered as reference and Φij is the component i, j of the transition matrix calculated by one of the methods. It gives a rough idea of the number of the common significant figures retained after the computations.The first method considers the pure Keplerian motion and the second one considers the Keplerian motion and theJ2 effect. A whole day of integration is divided in intervals of 1 to 60 seconds, which produces from 86,400 to 1,440 steps of integration, respectively, as shown in Table 1. The total processing time of each method is also shown in Table 1, although it depends on the code, the computer, and the programmer skills (the evaluation used a regular PC Pentium II processor, 512 MB, running FORTRAN codes). However, it illustrates roughly the comparative expected performance. One can note that the CPU time difference at any step of integration and for any method is very small. The first method is slightly faster than the second one for small steps of time, like 1 and 10 seconds.Table 1
Processing time for analytical calculation of the transition matrix.
Δt (s)StepsCPU TIME(s)First methodSecond method186,40019.826.2108,6402.12.7302,8800.70.9601,4400.40.4To analyze the accuracy of the methods, six comparisons are done, using two kinds of orbits: one circular and one elliptical. The transition matrix generated by each of the methods is compared with the reference generated by the RK78. The circular orbit is from the Topex/Poseidon satellite [21] and the elliptical one is from the Molniya satellite [12]. All those tests are performed for a period of 24 hours. The data of these orbits are shown in Table 2.Table 2
Keplerian parameters of the test orbits.
ParameterTopex/PoseidonMolniyaSemiaxis major7714423.46 m26563000.0 mEccentricity1.13458×10-40.75Inclination66°.03963°.435Longitude of node236°.720°.0Argument of perigee102°.83270°.0Mean anomaly153°.540°.0Besides, the RK4 (Runge-Kutta of fixed 4th order) numerical method for computing the transition matrix is also implemented and the transition matrix generated by the RK4 is also included in the comparison. The RK4 used the same dynamical model of the reference RK78. The results of those three comparisons for the circular orbit are shown in Table3, and the results for the elliptical orbit are shown in Table 4. The errors are shown in terms of per step mean and standard deviations considering the whole day sample. The comparison is done for only these two kinds of orbits, where the Keplerian approximation for the transition matrix provides already reasonable accuracy. For large time interval, the second method is always disadvantageous because it is a local numerical approximation, as depicted in Tables 3 and 4.Table 3
Comparison for a circular orbit (Topex).
Δt (s)Global errorRK4First methodSecond method13.7×10-3±6.0×10-31.2×10-6±1.3×10-56.3×10-8±8.9×10-7103.8×10-3±8.3×10-31.0×10-4±5.4×10-45.0×10-6±4.0×10-5303.9×10-3±9.5×10-37.6×10-4±2.1×10-33.5×10-5±1.4×10-4603.6×10-3±4.9×10-32.7×10-3±5.5×10-31.3×10-4±4.0×10-43003.2×10-3±1.7×10-36.6×10-2±1.1×10-13.8×10-3±8.2×10-36003.1×10-3±1.4×10-32.3×10-1±2.3×10-12.1×10-2±2.5×10-2Table 4
Comparison for an elliptical orbit (Molniya).
Δt (s)Global errorRK4First methodSecond method13.7×10-4±9.8×10-48.8×10-7±2.3×10-56.2×10-8±1.7×10-6103.6×10-4±8.2×10-45.7×10-5±1.1×10-34.0×10-6±7.6×10-6303.5×10-4±7.4×10-42.4×10-4±2.7×10-31.7×10-5±1.9×10-4603.4×10-4±6.9×10-45.1×10-4±3.1×10-33.6×10-5±2.3×10-43002.7×10-4±3.8×10-41.5×10-2±2.3×10-11.6×10-3±1.3×10-26003.5×10-4±6.7×10-44.3×10-2±1.9×10-14.6×10-3±1.9×10-2


## 6. Conclusions

The goal of this research was to compare and choose the most suited method to calculate the transition matrix used to propagate the covariance matrix of the position and velocity of the state estimator, (e.g., Kalman filtering or least squares), within the procedure of the artificial satellite orbit determination. The methods were evaluated according to accuracy, processing time, and handling complexity of the equations for two kinds of orbits: circular and elliptical. The processing time of the second method is around 30% larger than the first one; however, this difference is not considered enough to harm the computer load in the orbit determination tasks. The second method is more accurate for short intervals, asΔt=1 and 10 seconds, for the orbits considered. For other intervals of propagation, the first method shows to be more stable in the sense of keeping the same accuracy regardless of the step size. The equations of the second method are easier to be handled, which means that it is possible to include more perturbations easily. However, the first method is a closed analytical solution for the two-body problem only. The second method has no singularity and no restriction about the type of orbit; the first one is optimized for elliptical orbits (not optimized for parabolic and hyperbolic orbits). The second method performs an approximation in Taylor’s series, whereas in the first there is an analytical solution for the Keplerian motion. Therefore, for small intervals of time (1 to 10 seconds) and when one expects more accuracy, it is recommended to use the second method, since the CPU time does not overload excessively the computer in the orbit determination procedure. For larger intervals of time and when one expects more stability on the calculation, it is recommended to use the first method.



---

*Source: 768973-2011-10-29.xml*