1.

INTRODUCTION

Kepler’s equation (KE), which describes how a body moves under the influence of gravity, is derived in orbital mechanics. In his book Astronomia Nova [1], Johannes Kepler first discovered it. The elliptical KE is

where M, E and e designate mean anomaly, eccentric anomaly and eccentricity, respectively. Both M, E are fundamental parameters for determining the position of a moving celestial body in an elliptical orbit. KE is a transcendental equation because it involves a sine function. The exact analytical solution is unknown. It’s simple to calculate M for a given value of E. However, because there is no closed-form solution, the inverse issue, which involves finding E while M and e is known, can be far more difficult. Usually, E needs to be estimated by series expansions or numerical methods. Kepler himself approximate his equation by simple iteration in 1621 in his book Epitome of Copernican Astronomy [1]. KE is one of the core equations and has a lot of applications in orbital mechanics, therefore even though many academics have developed several ways to solve it, this subject continues to draw attention. For KE, finding a simple, accurate, and analytical solution is still of practical importance.

In this paper, the inverse problem is transformed to a machine learning (ML) problem, more exactly a supervised learning problem. By solving the supervised learning problem, a new solution is learned from the pre-calculated data. The proposed method, called ML-based method, takes the advantage of the great flexibility of neural networks (NNs). The approach is appropriate for extensive and quick orbit propagation since the complexity of the suggested algorithm is constant and independent of the eccentricity and transition time.

The rest of this paper is structured as follows. Existing approaches to resolving KE are carefully compiled and reviewed in Section 2. The primary idea behind the suggested machine learning-based solution is described in Section 3. And in Section 4, the effectiveness of the suggested strategy is verified using numerical simulation. Finally, Section 5 provides a summary of the result.

2.

METHODS FOR SOLVING KEPLER’S EQUATION

Many academics have looked into the inverse problem of KE. Finding approaches to arrive at the solution with high accuracy and minimal computational expense is the aim. Figure 1 illustrates how the various approaches to solving KE can be categorized. In this section, some very efficient and classical methods are introduced firstly as shown in Table 1. In Section 4, several of these techniques will be contrasted with our new technique in terms of their effectiveness and precision.

Figure 1.

Classification of ways to solve KE.

Table 1.

Ways to solve Kepler’s equation.

Kepler himself created the first method for resolving his equation, and Newton’s method came next. The idea of Newton’s method [2, 3] (or Newton-Raphson’s Method) is to approximate the KE by the first two terms in a Taylor series expansion. By extending the series to the n-term, a generalized Newton’s method is obtained. When n =3, the generated solution has an equivalent form to Halley’s method [4]. When n =4, the generated solution is same to Danby’s method [5]. Moreover, truncating the series to the first-order leads to a method which is identical to Kepler’s method [6]. Thus, Kepler’s method can also be seen as a variant of Newton’s method. In order to finding E while M and e is known, the initial value of E should be searched first. Then the repetition is continued until accuracy is satisfied. On the other hand, Newton’s technique and its variations may diverge if the original guess is not sufficiently close to the solution. According to Danby [5], as the degree of convergence increases, so do the initial value’s sensitivity and the risk of divergence. Although Newton’s approach and its variations perform well close to the solution, they lack a feature that would allow them to converge worldwide. The convergence of Conway’s method [7] and Mortari’s method [8] is guaranteed. However, these methods, like Newton’s method and its variants, has different iteration steps for different M and e which means the runtime is dependent on M and e. Except the classification result as shown in Table 1, the methods can also be divided into two categories: single-point method and sequential method. KE is solved using the single-point method, but this method does not benefit from the fact that KE has been computed at the previous neighborhood point. While sequential method [9-11] calculates the present value using the value from the previous moment, making full use of the previous calculation. However, this kind of methods are only applicable to orbit propagation problems with fixed time steps. In addition, the propagation of the initial error causes such algorithms to be sensitive to the error of initial step. Expansion method, including Lagrange expansion [2], Fourier-Bessel expansion [2] and Adomian Decomposition (AD) Method [12, 13], expands the solution of KE into a polynomial consists of N terms. This kind of algorithms are more suitable for the case of small eccentricity, because when the eccentricity is large, more terms need to be reserved in order to improve the accuracy, which leads to low calculation efficiency. The two-step method is a non-iterative method which divides the problem of solving KE into two steps. First, calculate an initial estimate by using a technology, such as Bezier Curve [14] and series expansion [15]. Then correct the initial value to a high accuracy using a generalized Newton’s method which is applied only once, rather than in a loop. Data-based method [16, 17] determines a start value according to the pre-computed data and then using iterative method to correct it.

For a number of issues in Celestial Mechanics, KE must be solved numerous times. Due to this, it is crucial for these applications to efficiently compute mean anomaly. For instance, in trajectory planning problems of the relative motion, the true anomaly should be calculated repeatedly which involves solving KE. The accuracy and speed of trajectory planning are significantly impacted by the effectiveness and stability of the eccentric anomaly calculation solution. Therefore, in this paper a method based on machine learning with low computation cost and constant complexity is presented.

3.

MACHINE LEARNING-BASED SOLUTION

This section presents an innovative method for solving KE. E is a function of M and e, defined as E = h(M, e). The function’s precise formulation cannot be expressed in closed form. Since 2π radians are swept off per period T, the eccentric anomaly may be calculated for any M in an infinite interval [0,+∞) by

where [‧] denotes the remainder and 0 ≤ [M / 2π] ≤ 2π. In addition, the function E = h(M,e) is centrosymmetric about E = M = π, namely

Therefore, the eccentric anomaly in the whole interval may be derived as long as the function E = h(M, e) is computed in the interval 0 ≤ M ≤ π.

In supervised learning, a machine learning job, an input is mapped to an output using examples of input-output pairings. With each training example including an input item and the desired output value, or “labelled training examples,” a function is created. A supervised learning algorithm produces an inferred function by looking at the training data. The inferred function can be used to approximate the unknown function E = h(M,e) in the inverse problem.

3.1

Feedforward neural network

Artificial NNs, or NNs for short, are a mathematical algorithmic model for distributed parallel information processing that replicates the behavioral properties of animal neural networks. Such networks rely on the system’s complexity to process information by altering the interaction between a large number of internally linked nodes.

Universal approximation theory (UAT) [18] pointed out that any bounded and regular function from one finite-dimension space to another can be approximated by an ordinary multilayer feedforward NN with any required accuracy. The network must have a sufficient number of neurons with a linear output layer and at least one hidden layer.

As shown in Figure 2, a feedforward NN may be used to estimate an unknown function h(M,e) and produce an analytical solution with high precision. For the approximation problem of the function E = h(M,e), there are two inputs and one output. One may think of the network as a composite function g(M,e). Thus, the estimate of eccentric anomaly is designed as

Figure 2.

The NN for solving KE.

In hidden layers, one of the squashing functions, the hyperbolic tangent sigmoid transfer function, is used in this research as the active function

and in the output layer, the activation function is

3.2

Generation of training data

Due to the lack of an accurate solution to the eccentric anomaly, we are unable to get the exact value E when given M and e. Calculating M when given E and e is, fortunately, simple. As a result, Algorithm 1 is made to acquire the training data.

3.3

Final adjustment

UAT shows that given enough neurons, we can obtain solutions to KE with arbitrary accuracy. Too many neurons can make it difficult to train the network. Therefore, in order to reduce the training difficulty of the network, a small network is chosen as the learning model in this paper. In this case, although an estimate very close to the true value E can be obtained, the accuracy may not be high enough. To improve the accuracy, a single-step adjustment algorithm is applied

The single-step adjustment algorithm designed by Halley’s method [4] is as follows

To increase the solution’s accuracy, a single-step adjustment is performed. The analytical solution’s complexity does not change because there is only one iteration.

4.

NUMERICAL SIMULATION

As of October 2020, one hundred and ninety-nine elliptical orbit objects with eccentricity greater than 2 were discovered in the publicly accessible data (www.space-track.org). Excentricity is 0.898411 at its highest level. Therefore, the greatest limit of eccentricity, according to our assumptions, is 0.9. Decompose E∈[0, π] and e∈ [0, 0.9] 2 into smaller intervals of equal size using eighty points respectively. 6400 training examples are obtained using Algorithm 1.

In this paper, we verify the performance of the proposed algorithm using a three-layer neural network that contains two hidden layers, each with nine neurons. Using Matlab toolbox, train the network with 6400 training examples and the residuals of the network at the examples are shown in Figure 3.

Figure 3.

Errors in computed E.

In practical engineering application, the inverse problems are solved with given eccentricities for different mean anomalies. Therefore, we take the max error in computed E under a given eccentricity, as shown in Figure 4, as index of accuracy of the algorithm. As shown in Figure 4 the algorithm proposed in this paper has high accuracy with final adjustment. When e = 0.89, the worst accuracy of 6.4×10^–11 is obtained. For e <0.7 the ML-based method has a precision of 10^–15.

Figure 4.

Maximum of errors in computed E under given eccentricities.

The series methods and the ML-based method can be classified into one category based on the fact that these two kinds of methods both use known models to approximate unknown functions. Figures 5 and 6 show that for small eccentricity the Lagrange method and Fourier-Bessel method can obtain same precision with ML-based method even with a few terms. For large eccentricity the ML-based method has a higher precision. Additionally, the Lagrange series diverges when it surpasses 0.662743419, which suggests that adding additional terms produces poorer outcomes. For large eccentricity the number of items of Lagrange method and Fourier-Bessel method is greatly increased. This makes Lagrange method and Fourier-Bessel method very computational expensive as shown in Table 2. Table 2 shows that ML-based method is more efficient than the Lagrange method and Fourier-Bessel method. With increasing of the number of terms used for improving solutions’ accuracy, the mean runtime of ML-based method is nearly constant and very small. However, the mean runtime of Lagrange method and Fourier-Bessel method increased dramatically.

Figure 5.

Comparison of the Lagrange method’s accuracy with the ML-based approach.

Figure 6.

Comparison of the Fourier-Bessel technique’s accuracy with the ML-based approach.

Table 2.

Mean runtime of Lagrange method, FB method and ML-based method.

5.

CONCLUSION

In this paper a non-iterative method based on machine learning for solving KE is presented. Based on neural networks’ tremendous flexibility, this technique can offer a high precision solution to KE. It can be seen as an analytical solving tool because the new method does not require any iterative computation or numerical integration. It also avoids the shortcoming of sensitivity of initial guess in some classical methods for solving KE. Numerical accuracy tests show that, the new algorithm behaves better than the most existed methods both in computational efficiency and accuracy. The approach outlined in this work, aside from the elliptical instance, may also be used to parabolic case and hyperbolic case.