Comparison of Numerical Orbit Integration between Runge-Kutta and Adams-Bashforth-Moulton using GLObal NAvigation Satellite System Broadcast Ephemeris

Abstract

Numerical integration is necessary for satellite orbit determination and its prediction. The numerical integration algorithm can be divided into single-step and multi-step method. There are lots of single-step and multi-step methods. However, the Runge-Kutta method in single-step and the Adams method in multi-step are generally used in global navigation satellite system (GNSS) satellite orbit. In this study, 4th and 8th order Runge-Kutta methods and various order of Adams-Bashforth-Moulton methods were used for GLObal NAvigation Satellite System (GLONASS) orbit integration using its broadcast ephemeris and these methods were compared with international GNSS service (IGS) final products for 7days. As a result, the RMSE of Runge-Kutta methods were 3.13m and 4th and 8th order Runge-Kutta results were very close and also 3rd to 9th order Adams-Bashforth-Moulton results. About result of computation time, this study showed that 4th order Runge-Kutta was the fastest. However, in case of 8th order Runge-Kutta, it was faster than 14th order Adams-Bashforth-Moulton but slower than 13th order Adams-Bashforth-Moulton in this study.

1. INTRODUCTION

A satellite navigation system consists of ground, space, and user segments. All satellite navigation systems have the same system configuration. The ground segment is then composed of elements required for monitoring and control for overall system operation. The ground segment performs the following functions: satellite status monitoring, maintenance plan establishment for system availability and minimization of performance degradation, satellite ephemeris and clock generation. Each system has a different number of facilities to perform these functions.

The ground segment utilizes system information such as satellite orbit status and maintenance plan to predict the satellite orbit and generates navigation message and transfers the message to the satellite, thereby making users to receive up-to-date information in real time. To predict the satellite orbit, a numerical integration algorithm is needed, by which satellite acceleration is predicted and the orbit can be predicted by applying the predicted acceleration to the velocity and position at the previous time.

In relation to the satellite orbit determination, Whalley (1990) developed a satellite orbit determination algorithm, compared the baseline analysis performance using global positioning system (GPS) broadcast ephemeris and precise ephemeris, and employed a 4th order Runge-Kutta and 9th order Adams predictor-corrector algorithm for the satellite orbit determination. In addition, Hein et al. (1997) employed the 8th order Runge-Kutta and 12th order Adams-Cowell algorithm for inclined geosynchronous orbit (IGSO) simulations of Europe regions. Su (2000) performed simulations using the 8th order Runge-Kutta and 9th order Adams-Cowell algorithm to analyze the strengths and drawbacks of geostationary earth orbit (GEO), ISGO, and medium earth orbit (MEO) determination. Bae (2009) modeled geopotential, third body, and solar radiation pressure for the acceleration of GPS satellites to analyze GPS satellite precise orbit determination and error characteristics and compared the results with the precise ephemeris provided by the international GNSS service (IGS), in which the 12th order Adams-Bashforth-Moulton algorithm was employed for acceleration integration.

In relation to the multi-step integration algorithm coefficient, Maury & Segal (1969) proposed various numerical integration algorithm coefficients up to the 16th order including Adams for round-off error of the algorithm and improvements on speed and accuracy, and Kirkpatrick (1976) proposed the Adams-Bashforth-Moulton coefficient up to the 20th order. However, the study explained that a higher order did not necessarily guarantee higher accuracy. Maury & Segal (1969) and Kirkpatrick (1976) denoted the coefficients of all orders as a number in their literature, whereas Seidu (2011) provided an equation by which Adams-Bashforth-Moulton coefficients could be calculated thereby being able to easily apply higher order coefficients.

In relation to the comparison between numerical integration algorithms, Fathoni & Wuryandari (2015) applied the calculation algorithms of acceleration, velocity, and position to Euler, Heun, 4th order Runge-Kutta, and 3rd order Adams-Bashforth-Moulton using gravity force, spring damper force, and aerodynamic force. Their calculation using C++ took about 0.055 seconds and the fastest computation speed was Euler followed by Heun, 4th order Runge-Kutta, and 3rd order Adams-Bashforth-Moulton.

Bate et al. (1971) and Xu (2008) reported that RungeKutta and Adams algorithms were typically used in numerical orbit integration and multi-step algorithm was faster than single-step algorithm, but the study by Fathoni & Wuryandari (2015) reported a different result. In addition, most previous studies applied numerical integration algorithms to GPS as a fixed order. Thus, the present study implemented various orders of Runge-Kutta and Adams-Bashforth-Moulton algorithms, and compared the accuracy and calculation time between numerical integration algorithms by applying the algorithms to GLObal NAvigation Satellite System (GLONASS) broadcast ephemeris.

2. GLONASS EPHEMERIS

The GLONASS broadcast ephemeris provides satellite position and velocity and perturbing accelerations due to the sun and the moon, by which a satellite position at a preferred time can be calculated through numerical integration (GLONASS ICD 2008). The calculation method of GLONASS satellite position has three equations, and the RungeKutta algorithm (GLONASS ICD 2016) can be used for the calculation. Table 1 presents the numerical integration errors of three satellite position calculation methods provided in GLONASS ICD (2016).

Table 1. Position prediction errors at various intervals (m) (GLONASS ICD 2016).

Algorithm example	Integration interval
	5 min	10 min	15 min	4 h
Precise Simplified Long-term	0.13 0.03-0.42 0.03-0.42	0.18 0.04-0.56 0.04-0.56	0.25 0.05-0.77 0.05-0.77	>30 >100 0.25-1

 

As presented in Table 1, the simplified and long-term methods have the same error except for four-hour integration interval. Since GLONASS broadcast ephemeris is given at every 30 min and integrated in every 15 min from the reference time of ephemeris in general, this study employed a simplified method. Eq. (1) presents a simplified equation to calculate GLONASS satellite acceleration.

\(\begin{aligned} \frac{d V_{x}}{d t} &=-\frac{G M}{r^{3}} x-\frac{3}{2} J_{2}^{0} \frac{G M \cdot a_{e}^{2}}{r^{5}} \times\left(1-5 \frac{z^{2}}{r^{2}}\right)+\omega_{E}^{2} x+2 \omega_{E} V_{y}+\ddot{x} \\ \frac{d V_{y}}{d t} &=-\frac{G M}{r^{3}} y-\frac{3}{2} J_{2}^{0} \frac{G M \cdot a_{e}^{2}}{r^{5}} y\left(1-5 \frac{z^{2}}{r^{2}}\right)+\omega_{E}^{2} y-2 \omega_{E} V_{x}+\ddot{y} \\ \frac{d V_{z}}{d t} &=-\frac{G M}{r^{3}} z-\frac{3}{2} J_{2}^{0} \frac{G M \cdot a_{e}^{2}}{r^{5}} \mathrm{z}\left(3-5 \frac{z^{2}}{r^{2}}\right)+\ddot{z} \end{aligned}\)       (1)

Here, GM is the geocentric gravitational constant, which is (398600441.8 ± 0.8) × 10⁶ m³/s², r refers to the distance between Earth center and satellite, x,y,z refer to the satellite position, \(j^0_2\) refers to the second degree zonal coefficient of normal potential, which is 1082625.75 × 10^-9, a_e refers to semi-major axis of the Parametry Zemli 1900 (PZ-90) Earth’s ellipsoid, which is 6378136 m, ω_E refers to the mean angular velocity of the Earth relative to the vernal equinox, which is 7.2921151467 × 10^-5 rad/s, V_x, V_y refer to the satellite velocity, and \(\ddot{x}\), ÿ, \( \ddot{z}\) refer to the perturbing accelerations of the sun and the moon provided by the broadcast ephemeris(GLONASS ICD 2008, 2016).

3. NUMERICAL INTEGRATION ALGORITHM

3.1 Runge-Kutta Method

Runge-Kutta is a typical single-step numerical integration method, by which orbit integration can be done using a single initial value. Runge-Kutta is relatively simple and easy to be implemented, and an integration step-size can be easily changed. However, it is difficult to find an appropriate integration step-size because there is no simple way to determine the truncation error (Bate et al. 1971).

Xu (2008) derived an equation through assumptions because it could not have a single solution as the number of unknown values was thirteen and the number of equations was seven in the expansion of the 4th order Runge-Kutta equation using the Taylor series. Accordingly, it was verified that despite of the same 4th order Runge-Kutta algorithm, it was different from the algorithm proposed by Bate et al. (1971). In this study, the 4th order Runge-Kutta algorithm proposed by Bate et al. (1971) was used, which is presented in Eq. (2). Some of the algorithm was modified in Eq. (2) to explain the equation conveniently according to this study.

\(\begin{aligned} a_{\text {integration }} &=\frac{1}{6}\left(K_{1}+2 K_{2}+2 K_{3}+K_{4}\right) \\ K_{1} &=F\left(P_{0}, V_{0}, a_{3 r d}\right) \\ K_{2} &=F\left(P_{K_{1}}, V_{K_{1}}, a_{3 r d}\right) \\ K_{3} &=F\left(P_{K_{2}}, V_{K_{2}}, a_{3 r d}\right) \\ K_{4} &=F\left(P_{K_{3}}, V_{K_{3}}, a_{3 r d}\right) \end{aligned}\)       (2)

Here, a_integration refers to the integrated acceleration value, F refers to Eq. (1), P₀, V₀ refer to the initial position and velocity used in the integration, and a_3rd refers to the perturbing accelerations due to the sun and the moon. In addition, P_Kn, V_Kn refer to the position and velocity values calculated using K_n, which follow the position, velocity, and acceleration equations as presented in Eq. (3).

\(\begin{array}{l} P_{K_{n+1}}=P_{0}+V_{0} K_{n} t+\frac{1}{2} K_{n} \times \alpha(\beta t)^{2} \\ V_{K_{n+1}}=V_{0}+K_{n} \beta t \end{array}\)       (3)

Here, t refers to the integration step-size, which is a second unit of time in this study, and α, β are presented in Table 2.

Table 2. 4th Runge-Kutta coefficient of Eq. (3).

For	\(\alpha \)	\(\beta \)
\( K_2 \) \( K_3 \) \( K_4 \)	1 1 1	1/2 1/2 1

 

t-integrated position and velocity can be calculated using Eq. (3), K_n refers to a_integration, and α, β refer to 6 and 1/6, respectively.

The 8th order Runge-Kutta was the same equation proposed by Su (2000) and Xu (2008), but it was expressed in different manner. The equation proposed by Su (2000) is presented in Eq. (4), and some of the equation was omitted.

\(\begin{array}{l} y_{n+1}=y_{n}+\frac{h}{840}\left(41 K_{1}+27 K_{4}+272 K_{5}+27 K_{6}+216 K_{7}+216 K_{9}+41 K_{10}\right) \\ K_{1}=F\left(t_{n}, y_{n}\right) \\ K_{2}=F\left(t_{n}+\frac{4}{27} h, y_{n}+\frac{4 h}{27} K_{1}\right) \\ K_{3}=F\left(t_{n}+\frac{2}{9} h, y_{n}+\frac{h}{18}\left(K_{1}+3 K_{2}\right)\right) \\ \vdots \\ K_{10}=F\left(t_{n}+h, y_{n}+\frac{h}{820}\left(1481 K_{1}-81 K_{3}+7104 K_{4}\right.\right. \\ \left.\left.\quad-3376 K_{5}+72 K_{6}-5040 K_{7}-60 K_{8}+720 K_{9}\right)\right) \end{array}\)       (4)

α, β in Eq. (3)m which calculate the position and velocity using the 8th order Runge-Kutta are presented in Table 3, and h-interated position and velocity are: in Eq. (3), K_n is y_n+1, and α, βare 840 and 1/840, respectively.

Table 3. 8th Runge-Kutta coefficient of Eq. (3).

For	\(\alpha \)	\(\beta \)
\( K_2 \) \( K_3 \) \( K_4 \) \( K_5 \) \( K_6 \) \( K_7 \) \( K_8 \) \( K_9 \) \( K_{10} \)	1 4 4 4 36 720 20 480 820	4/27 1/18 1/12 1/8 1/54 1/4320 1/20 1/288 1/820

 

As presented in Eqs. (2, 4), Runge-kutta requires a number of K calculations. Thus, its computation speed is relatively slower than that of the multi0step because self-starting is possible (Bate et al. 1971, Su 2000, Xu 2008).

3.2 Adams-Bashforth-Moulton Method

The Adams algorithm, which is one of the multi-step integration methods, is divided into predictor and corrector. The predictor algorithm is called Adams-Bashforth and the corrector algorithm is called Adams-Moulton. Since AdamsMoulton can be calculated iteratively, an integration stepsize can be adjusted and the number of orders can also be adjusted when the truncation error is larger than the required error. In addition, since only one Eq. (1) calculation is done in the integration phase, its calculation speed is faster than single-step method (Bate et al. 1971, Xu 2008).

According to Wikipedia (2019), the first order AdamsBashforth is called Euler method, and the 1st-order and 2nd order Adams-Moulton algorithms are called backward Euler method and trapezoidal rule, respectively. Thus, AdamsBashforth-Moulton can be used from the 3rd order. Three initial values are needed to employ the 3rd order AdamsBashforth-Moulton algorithm, which is calculated using Runge-Kutta in general (Bate et al. 1971, Xu 2008).

The predictor algorithm is presented in Eq. (5), in which a_p refers to the predicted GLONASS acceleration, h refers to the integration step-size, which is a second unit of time in this study, n refers to the order number, and C_B refers to the coefficient for integration, and every order has different coefficient. a_i refers to the previous acceleration value for prediction, which requires n values. In addition, ai refers to the value calculated through Eq. (1). Note that it is not a_integration in Eq. (2) and y_n+1 in Eq. (4).

\(a_{p}=\mathrm{h} \sum_{i=1}^{n} C_{B} a_{i}\)       (5)

The corrector algorithm is presented in Eq. (6), in which a_c refers to the corrected GLONASS acceleration and C_M refers to the coefficient for integration, which is different from C_B. a_i refers to the acceleration value for correction, which requires a predicted value using Adams-Bashforth and n-1 previous values.

\(a_{c}=h \sum_{i=1}^{n} C_{M} a_{i}\)       (6)

The position and velocity can be calculated by using a_p or a_c and applying it to Eq. (3), in which α=1, β=1. In addition, the truncation error can be calculated using the difference between a_p and a_c, and the calculation was iterated in this study when the difference of the satellite position is larger than 1×10^-6 m. The number of iterations in most orders is two in maximum.

4. COMPARISION RESULTS

In thin section, GLONASS satellite positions were calculated using the implemented algorithm, and the performance was compared with the IGS final products. Kouba (2015) reported that the satellite coordinates were provided refer to the center of mass in the IGS products but the broadcast ephemeris like the GPS was provided refer th the phase center, which required correction. However, according to GLONASS ICD (2016), GLONASS broadcast cphemeris is provided refer to the center of mass. Thus, the coordinate correction process can be omitted. Because the IGS precise ephemeris provided the satellite position in every 15 min based on global positioning system time (GPST), it was re-calculated in every second using the 9th order Lagrange interpolation of GLONASS time. The analysis period was from 2019 day of year (DOY) 195 to 201. Since Nos. 6 and 25 satellites were not present in the products, and No. 26’s health flag was 1 (unhealthy) in the broadcast ephemeris, they were excluded from the analysis.

4.1 4th and 8th Order Runge-Kutta Method

The GLONASS broadcast ephemeris provides the parameter in every 30 min and performs forward and backward computations of 15 min as shown in Fig. 1. Here, there will be an overlapped portion as shown in 30 min in Fig. 1.

Fig. 1. GLONASS broadcast ephemeris integration general scheme.

The 30-min portion in Fig. 1 was calculated using the forward and backward manner, and their accuracy was compared. To do this, the integration was performed using the 4th and 8th order Runge-Kutta. There are two methods in the integration using Runge-Kutta as shown in Fig. 2. In relation to this, Habrich (1999) reported that significant errors occurred as the integration step-size was larger, but the difference between one and ten seconds was similar at a level of several tens of centimeters. Thus, this study employed one and five seconds step-size that was used in CASE 1 of Fig. 2 method to perform the integration and compared the results.

Fig. 2. Two different integration method of GLONASS broadcast ephemeris.

Table 4 presents the comparison results with the IGS final products by the root mean square error (RMSE). The maximum, minimum, mean, and standard deviation of 23 satellites used in the analysis are presented. Although 48 ephemeris were provided daily as the comparison data for one week, DOY 195 0:0:0 in the backward and DOY 2020:0:0 in the forward were excluded. Thus, 335 ephemeris per satellite were used. As presented in Table 4, the forward method exhibited slightly better results than that of the backward method, and the 4th order results were slightly better than those of the 8th order. Excluding the maximum value of the forward, the results of one second step-size integration were better than those of the five seconds stepsize integration at a level of centimeter.

Table 4. Forward and backward calculation RMES result of 4th and 8th Runge-Kutta method compared with IGS final product.

	4th Runge-Kutta				8th Runge-Kutta
Scheme	Forward		Backward		Forward		Backward
Step-size	1s	5s	1s	5s	1s	5s	1s	5s
Max (m) Min (m) Ave. (m) Std. (m)	6.51 2.29 3.17 1.06	6.41 2.31 3.17 1.04	6.95 2.41 3.21 1.13	7.07 2.40 3.22 1.14	6.51 2.29 3.17 1.06	6.42 2.33 3.18 1.04	6.95 2.41 3.22 1.13	7.08 2.42 3.24 1.14

 

The weekly RMSE analysis results after integration in every second showed that X, Y, Z, and 3D were 1.80 m, 1.79 m, 1.82 m, and 3.13 m on average, respectively, and no significant deviation was found among X, Y, and Z. The results of the 4th and 8th order Runge-Kutta were different at a level of millimeter. In relation to the comparison between broadcast and precise ephemeris of GLONASS, Lee (2012) performed 24-hour analysis in every second and the analysis results showed that X, Y, and Z were 5.11 m, 5.09 m, and 4.13 m, respectively. In comparison to this study result, the accuracy of GLONASS broadcast ephemeris has significantly improved. Fig. 3 shows the RMSE results for each satellite of the 8th order Runge-Kutta, in which the X-axis and Y-axis represent the satellite number and m-unit RMSE, respectively.

Fig. 3. 8^th Runge-Kutta integration RMSE compared with IGS final product for 7 days.

4.2 Various Orders of Adams-Bashforth-Moulton Method

By using the calculation equation of Adams-BashforthMoulton coefficient proposed by Seidu (2011), AdamsBashforth and Adams-Moulton can calculate up to 81th order and 82th order, respectively. In this study, up to the 81th order was applied. Fig. 4 shows the 3D RMSE results, which are compared with the IGS final products by integrating in every second. The 8th order Runge-Kutta was used to calculate the initial value, and the initial value was excluded in the RMSE calculation.

Fig. 4. Adams-Bashforth-Moulton integration RMSE compared with IGS final product for 7 days.

As shown in Fig. 4, the 3D RMSE was around 3.13 m from the 3rd to 9th orders, which was very similar to that of the 4th and 8th order Runge-Kutta, while the 3D RMSE at the 10th order was 3.14 m, which was slightly different. The error increased gradually from the 10th order, and the 3D RMSE was the largest as 11 m at the 13th order. In relation to the error increase, the condition number started to increase from the 11th order when using the computation equation of AdamsBashforth-Moulton coefficient proposed by Seidu (2011). The condition number increased as the order number increased, and the condition number at the 14th order coefficient was larger than that of the 13th coefficient. However, the RMSE in the 14th order was smaller than that of the 13th order, which suggests that additional analysis on other factors are necessary in addition to the condition number.

4.3 Computation Time

According to Bate et al. (1971), in general, the single-step integration is slower than the multi-step integration. The 4th order Runge-Kutta requires at least five times of computation processes for integration whereas the 8th order requires eleven times. For the 3rd order of Adams-Bashforth-Moulton, at least six times of computations are required, and iterations may occur depending on truncation error condition and it will take time of the previous acceleration data search. Thus, it is considered slower than the 4th order Runge-Kutta. To compare the computation time of the Runge-Kutta and Adams-Bashforth-Moulton algorithms, an integration time of one satellite for 30 min was computed, which is shown in Fig. 5. MATLAB was used to integrate every second in the computation. For Adams-Bashforth-Moulton, initial value computation and coefficient computation times were also included. Considering the degradation of central processing unit utilization rate assigned to MATLAB due to background programs running intermittently in the Windows Operating Systems, a mean value was used after measuring 10 times. The X-axis in Fig. 5 represents the order number of AdamsBashforth-Moulton, and the blue and purple lines represent the 4th order and 8th order Runge-Kutta, respectively.

Fig. 5. Computation time of 4^th and 8^th Runge-Kutta and various order of Adams-Bashforth-Moulton.

As shown in Fig. 5, the computation time of the 4th order Runge-Kutta was 43 milliseconds (ms), which was the fastest, and the computation time of the 8th order Runge-Kutta was 166 ms, which was almost three times slower than that of the 4th order Runge-Kutta. The computation time of the 3rd order Adams-Bashforth-Moulton was 80 ms. These figures verified that as the order number increased, the computation time increased by 3 ms - 6 ms. However, the computation times at the 9th and 11th orders were 28 ms and 17 ms, respectively, which exhibited somewhat longer time than those of other order numbers. The computation times at the 13th and 14th order Adams-Bashforth-Moulton were 162 ms and 168 ms, respectively, which revealed that the computation time of the 8th order Runge-Kutta was slower than that of the 13th order Adams-Bashforth-Moulton.

The difference between the 4th and 8th orders of RungeKutta, which was used to calculate the initial value of AdamsBashforth-Moulton, was 2 ms on average, which was very small, and the computation time at the 68th order was the largest as 15 ms.

5. CONCLUSIONS

This study implemented the Runge-Kutta and Adams-Bashforth-Moulton algorithms, which are the widely used numerical integration algorithms in single-step and multistep. The algorithms were applied to GLONASS broadcast ephemeris, and the results were compared with IGS final products. The accuracy of the two algorithms had no significant difference as their 7-day RMSE mean was 3.13 m. However, the accuracy of Adams-Bashforth-Moulton tended to degrade from the 10th order. The computation time of the 4th order Runge-Kutta was about twice faster than that of the 3rd order Adams-Bashforth-Moulton. Bate et al. (1971) reported that the single-step was slower than the multi-step in general. However, such conclusion was based on the result applied to predictor or corrector algorithm only. Accordingly, when performing the 2nd order Adams-Bashforth only, its computation time was 37 ms, which was faster than that of the 4th order Runge-Kutta. The accuracy was poor than that of the 4th Runge-Kutta, but the difference was very small only at a level of mm. However, an algorithm such as Runge-Kutta that calculates an initial value is needed to employ the Adams algorithm, and the computation time is similar with 4th order Runge-Kutta. Thus, this study verified that using the 4th order Runge-Kutta was considered sufficient for GLONASS broadcast ephemeris integration. This study completed the verification of the algorithm implemented using GLONASS broadcast ephemeris. The study results are expected to be utilized in the modeling of accelerations due to Earth’s gravitational field, celestial gravity, and solar radiation pressure, as performed by Bae (2009), and used as a basis algorithm for orbit prediction through the integration of GNSS satellites.

AUTHOR CONTRIBUTIONS

Conceptualization, D. W. Lim and S. Chun; methodology, D. W. Lim and E. Son; software, E. Son and J. Ahn; validation, E. Son, J. Ahn and M. Shin; formal analysis, E. Son, D. W. Lim and M. Shin; writing—original draft preparation, E. Son, D. W. Lim and S. Chun.; writing—review and editing, D. W. Lim, J. Ahn and M. Shin; visualization, E. Son

CONFLICTS OF INTEREST

The authors declare no conflict of interest.