# Autonomous Orbit Determination for Lagrangian Navigation Satellite Based on Neural Network Based State Observer.

1. Introduction

First the dynamical model of ERTBP with perturbation is introduced. Then we design a neural network based state observer to determine the orbit of Lagrangian satellites. Afterwards the stability of the observer is proved. Finally four Lagrangian satellites are chosen to validate the effectiveness of this AOD method.

2. Elliptic Restricted Three-Body Problem with Perturbation

As shown in Figure 1, [P.sub.1] and [P.sub.2] are the primaries in the three-body system. [P.sub.1] and [P.sub.2] are in elliptical orbits. P is the third body which is vanishingly small compared to the two primaries. Similar to the CRTBP, we study the motion of P in a rotating coordinates.

Let C-XYZ represent the barycentric synodic coordinate depicted in Figure 1; the x-axis of this frame is along the radius vector, which connects the primaries, positive in the direction pointing from [P.sub.1] to [P.sub.2]. The y-axis of this frame is perpendicular to the x-axis, positive in the direction of the motion of [P.sub.2]. The XY frame rotates with an angular velocity equal to the instant motion of [P.sub.2] with respect to [P.sub.1]. The z-axis is perpendicular to the orbital plane of the primaries.

The equations of motion in ERTBP are dimensionless. The dimensionless units are as follows:

[M] = [m.sub.1] + [m.sub.2]

[L] = [a (1 - [e.sup.2])]/[1 + e cos f]

[T] = [[[L.sub.3]/[G ([m.sub.1] + [m.sub.2])]].sup.1/2] = [square root of (1 + e cos f)]/[??], (1)

where [m.sub.1] and [m.sub.2] are the masses of the two primaries. a, e, respectively, refer to the semimajor axis and eccentricity of the two primaries' elliptical relative revolving orbit. f is the true anomaly of the secondary on the elliptic orbit. Then the motion of the Lagrangian satellite in the barycentric synodic coordinate system is governed by the following [18]:

X" - 2Y' = 1/[1 + e cos f] [partial derivative][OMEGA]/[partial derivative]X

Y" + 2X' = 1/[1 + e cos f] [partial derivative][OMEGA]/[partial derivative]Y

Z" = 1/[1 + e cos f] [partial derivative][OMEGA]/[partial derivative]Z. (2)

f is taken as the time-like independent variable. The first and second derivatives of the coordinate with respect to f are calculated as

X' = dX/df,

X" = [d.sup.2]X/d[f.sup.2]. (3)

[OMEGA] is the pseudo-potential function of the three-body problem; it is described as follows:

[OMEGA] = 1/2 [[X.sup.2] + [Y.sup.2] + [Z.sup.2] + [mu] (1 - [mu])] + [1 - [mu]]/[r.sup.1] + [mu]/[r.sub.2], (4)

where [mu] = [m.sub.2]/([m.sub.1] + [m.sub.2]), [r.sub.1] = [square root of ([(X + [mu]).sup.2] + [Y.sup.2] + [Z.sup.2])], [r.sup.2] = [square root of ([(X - 1 + [mu]).sup.2] + [Y.sup.2] + [Z.sup.2])].

In this study, only the motions around the collinear libration point [L.sub.1] or [L.sub.2] are investigated. Therefore the origin of the coordinate can be moved from the barycenter of system to the interested libration point for convenient analysis. As shown in Figure 2, the instantaneous distance between the libration point and its closest primary which is denoted as [gamma] is chosen as the new length unit. The new reference frame is defined as the [L.sub.1]- or [L.sub.2]-centered synodic reference frame. The directions of the axes of this new coordinate system are coincident with the barycentric synodic reference frame.

Let (x, y, z, x', y', z') represent the state variables in [L.sub.1]- or [L.sub.2]-centered synodic reference frame. The transformation of coordinates between the barycentric synodic frame and the [L.sub.1]- or [L.sub.2] -centered synodic system is represented as follows [19]:

X = [gamma] (x [- or +] 1) + 1 - [mu]

Y = [gamma]y

Z = [gamma]z, (5)

where the upper sign refers to the [L.sub.1] case and the lower one refers to the [L.sub.2] case. The linearized equations of motion in the [L.sub.1]- or [L.sub.2]-centered synodic system can be formulated as follows [20]:

x" - 2y' - (1 + 2[c.sub.2])x = [summation over (i[greater than or equal to]1)] [[(-e).sup.i] [cos.sup.i] f (1 + 2[c.sub.2]) x]

y" + 2x' - (1 - [c.sub.2])y = [summation over (i[greater than or equal to]1)] [[(-e).sup.i] [cos.sup.i] f (1 - [c.sub.2]) y]

z" + [c.sub.2]z = [summation over (i[greater than or equal to]1)] [[(-e).sup.i] [cos.sup.i] f (1 - [c.sub.2]) z], (6)

where [c.sub.2]([mu]) = (1/[[gamma].sup.3.sub.j])([mu] + (1 - [mu])[[gamma].sup.2.sub.j]/ [(1 [- or +] [[gamma].sub.j]).sup.3]) and [[gamma].sub.j] (j = 1, 2) is the instantaneous distance between Lj and its closest primary.

We define a new variable [bar.X] = [[x, y, z, x', y', z'].sup.T]; then (6) can be written as

X' = A[bar.X], (7)

where

[mathematical expression not reproducible] (8)

In addition to the gravitational force from the two primaries, several perturbations can affect the motion of Lagrangian satellite. When these perturbations are considered, the motion equation of the Lagrangian satellite can be described as

[bar.X]' = A[bar.X] + Bg ([bar.X]), (9)

where

[mathematical expression not reproducible] (10)

is the nonlinear perturbing force on per unit mass. And

[mathematical expression not reproducible]. (11)

3. The Observation

In this paper, satellite-to-satellite range is the observation. The relationship between the observation and the state is

[bar.Y] = [rho]([bar.X], t) = [square root of ([(x - [x.sub.2]).sup.2] + [(y - [y.sub.2]).sup.2] + [(z - [z.sub.2]).sup.2])], (12)

where [[x y z].sup.T] and [[[x.sub.2] [y.sub.2] [z.sub.2]].sup.T] are the position variable of Lagrangian satellite 1 and Lagrangian satellite 2, respectively, in [L.sub.1]- or [L.sub.2]-centered synodic reference frame.

Expanding (12) in a Taylor series with respect to estimated state yields

[mathematical expression not reproducible]. (13)

We get the linearized relation between observation and states by neglecting the higher order terms.

[mathematical expression not reproducible]. (14)

Now define the state estimate error as

[??] = [bar.X] - [??]. (15)

Further, the residual error of observation is defined as

[??] = [bar.Y] - [??]. (16)

Then we get

[??] = C[??], (17)

where

C = [[[??] - [x.sub.2]]/[rho] [[??] - [y.sub.2]]/[rho] [[??] - [z.sub.2]]/[rho] 0 0 0 ]. (18)

Usually we only can get the estimated value of the position of Lagrangian satellite 2; thus the matrix C can be calculated by using the following:

C = [[[??] - [[??].sub.2]]/[rho] [[??] - [[??[.sub.2]]/[rho] [[??] - [[??].sub.2]]/[rho] 0 0 0 ]. (19)

4. Design of the Adaptive Observer Based on Neural Network

The perturbation g([bar.X]) in (9) can be estimated using a neural network over a compact set [D.sub.[bar.X]] as follows [17]:

g ([bar.X]) = [W.sup.T][phi] ([bar.X]) + [epsilon] ([bar.X]) [parallel][epsilon][parallel] [less than or equal to] [epsilon]*, [bar.X] [member of] [D.sub.[bar.X]], (20)

where W is a matrix of unknown neural network weights, [phi]([bar.X]) is a known vector of bounded basis functions, and [epsilon]* is a uniform bound on the approximation error.

In order to estimate the state of Lagrangian satellite using the satellite-to-satellite range, an observer is designed as

[??] = A[??] + B [[??] ([??]) - [upsilon](f)] + K[??], (21)

where [upsilon](f) is a robust term that will be determined later. f is the time-like independent variable. K is a user-specified gain matrix. [??]([??])) is the estimated uncertainty vector which is given by

[??] ([??]) = [[??].sub.T][phi] ([??]). (22)

A weight estimate update [??] is designed as

[??] = F ([phi] ([??]]X) [[??].sup.T] - [sigma][??]), (23)

where F is a positive-definite symmetric matrices that will be determined later. And

[upsilon] (f) = -D [??]/[parallel][??][parallel] - [[epsilon].sub.max][??]. (24)

[phi]([??]) is calculated as follows [17]:

[phi] ([??]) = [1 [??]/[R.sub.[direct sum]] [??]/[R.sub.[direct sum]] [??]/[R.sub.[direct sum]] [??][??]/[R.sub.[direct sum]] [??][??]/[R.sub.[direct sum]] [??][??]/[R.sub.[direct sum]]]. (25)

Now we have completed the design of the observer. However when this method is used to determine the orbits of the Lagrangian satellite in a navigation constellation, the number of the Lagrangian navigation satellites must satisfy a constraint condition. In the following, we will derive the minimum number of the Lagrangian navigation satellites which is needed in this method.

In order to prove that the designed observer in (21) can estimate the state of Lagrangian satellites accurately, we define the errors in the neural network weights and the basis vector [phi] as follows:

[??] = W - [??]

[??] = [phi]([bar.X]) - [phi] ([??]). (26)

And a candidate Lyapunov functions is chosen as [17]

V = 1/2 [[??].sup.T]P[??] + 1/2 tr [[[??].sup.T][F.sup.-1]W], (27)

where P and F are positive-definite symmetric matrices. By differentiating (27), the Lyapunov function's derivative is illustrated as

[mathematical expression not reproducible]. (28)

We define

[A.sub.cl] = A - KC. (29)

Then it yields

[mathematical expression not reproducible]. (30)

If PB = [C.sup.T], it will yield

[mathematical expression not reproducible]. (31)

Here the Meyer-Kalman-Yakubovich (MKY) lemma [17, 21] is used to derive PB = [C.sup.T]. Noting that P ~ 6 x 6, B ~ 6 x 3, if PB = [C.sup.T], C should satisfy C ~ 6 x 3. It means that three satellite-to-satellite ranges should be provided. In other words, the navigation constellation should include four Lagrangian satellites which can generate three satellite-to-satellite ranges.

If four Lagrangian satellites construct the navigation constellation, C can be derived as [mathematical expression not reproducible], (32)

where [[[[??].sub.3] [[??].sub.3] [[??].sub.3]].sup.T] and [[[[??].sub.4] [[??].sub.4] [[??].sub.4]].sup.T] are the estimated position variable of Lagrangian satellite 3 and Lagrangian satellite 4, respectively.

Define matrix Q as

-Q = P[A.sub.cl] + [A.sup.T.sub.cl]P, (Q > 0). (33)

The Lyapunov function's derivative becomes

[mathematical expression not reproducible]. (34)

Since [??] = -[??] and [??] = F([phi]([??])[[??].sup.T] - [mathematical expression not reproducible]. (35)

From [17], it is noted that

tr [[W.sup.T][??][[??].sup.T]] [less than or equal to] [[??].sup.T]] [less than or equal to] [alpha] [parallel][??][parallel]. (36)

We also can get

[mathematical expression not reproducible]. (37)

[[parallel]x[parallel].sub.F] is calculated as the Frobenius norm.

Substituting (36)-(37) into (35), we can get

[mathematical expression not reproducible]. (38)

Thus, from (38), [??] < 0 when

[mathematical expression not reproducible], (39)

or

[mathematical expression not reproducible]. (40)

Note that parameters D, [[epsilon].sub.max], [sigma] are user selected so they can be appropriately chosen to minimize the estimation error.

5. Simulation and Analysis

In order to illustrate the effectiveness of the method introduced in this paper, we design a constellation around the Earth-Moon L1 libration point which includes four Lagrangian satellites. The initial states of the four satellites are listed in Table 1. The data are normalized. Quasiperiodic orbits are sensitive to initial states. Thus the initial states must be accurate enough. This is why the data in Table 1 retain many valid digits after the dot. We calculate quasiperiodic orbits with different initial states. The conditions in Table 1 may be different from the actual task situation. Therefore the initial conditions and measurement noise here are assumed conditions based on the characteristics of CRTBP.

Four quasiperiodic orbits are shown in Figure 3. The ephemeris of the Moon and Earth is obtained by DE405. In order to illustrate the AOD error, the "true trajectory" of Lagrangian satellites should be provided. In this paper, the "true trajectory" is calculated by using the precise dynamical equations as follows:

[??] = -[[mu].sub.3]/[R.sup.3] R - [n.summation over (i=1)] [[mu].sub.i] ([[DELTA].sub.i]/[[DELTA].sup.3.sub.i] + [R.sub.i]/[R.sup.3.sub.i]), (41)

where R is the position vector of the Lagrangian satellite in the J2000 coordinate. [[DELTA].sub.i] is the position vector of the ith celestial with respect to the Moon. [R.sub.i] is the position vector of the ith celestial with respect to the Earth. [[mu].sub.i] is the gravitational parameter of ith celestial. The quasiperiodic orbits are calculated by using multiple shooting method. The initial targets are obtained by using approximate analytic solutions. The initial targets are expressed in the barycentric synodic coordinate. Therefore, we should translate the initial targets into J2000 coordinate. Then the initial target is modified according to the task constraint, and the ideal quasiperiodic trajectory is obtained. Since the Lagrangian satellites in the Earth-Moon system is far away from the Earth, the perturbation caused by atmospheric resistance, temporal change in gravity, oblateness, and mass distribution is not considered. The main perturbations come from the gravity of the Sun and other planets. Therefore, the above equation is appropriate to describe the precise motion of the Lagrangian satellites.

The various tuning constants used in this scenario were chosen as follows:

[mathematical expression not reproducible]. (42)

The parameters are designed as [[epsilon].sub.max] = 1 x [10.sup.-12], D = 1x [10.sup.-9]. The step of AOD is 1 hour. Three scenarios with different initial state errors and observation errors as shown in Table 2 are analyzed. The gravitational attraction from the Sun is the only perturbation considered in the simulation.

For scenario 1, we give the simulation results of two Lagrangian satellites S1 and S2. The AOD errors in 30 days are presented in Figures 4(a) and 4(b). It can be seen that the AOD error of the both Lagrangian satellites is less than 10 m. The RMS of the AOD error of S1 in x-, y-, and z-axis are 2.3 m, 2.6 m, and 2.1 m, respectively. From Figures 5(a) and 5(b), we can see that this method can estimate the perturbation with small errors. The red curves are the theoretical calculating value of the Sun gravitational perturbation. The blue curves are the estimation of the Sun's gravity. The RMS of the estimate error of perturbation of S1 in x-, y-, and z-axis are 6.1 x [10.sup.-6] m/[s.sup.2], 5.9 x [10.sup.-6] m/[s.sup.2], and 4.3 x [10.sup.-7] m/[s.sup.2].

In Figures 6(a), 6(b), 7(a), and 7(b), we give the results for Scenario 2 and Scenario 3. As we can see, the maximum AOD errors increase with the initial state errors. All the stable AOD errors are less than 10 m. As shown in Figure 7(a), the AOD errors can converge to less than 10 m in 5 days, even the initial errors 100 m.

6. Conclusions

We analyzed the AOD for Lagrangian navigation constellation based on accurate dynamical mode. A neural network based method is applied to estimate the solar gravitational perturbation and the state of the Lagrangian satellites. The update rate of neural network weights is obtained by convergence of Lyapunov function. Therefore, the stability of the method introduced in this paper is proven. We derived that the constellation must include more than four satellites when this method is used. The simulation results show that this method can achieve a stable AOD error about 10 m only using satellite-to-satellite range as observation. Since the CRTBP is an approximate model for the Earth-Moon system, the AOD based on perturbed ERTBP is much useful for the future Lagrangian navigation constellation. We only test a Lagrangian navigation constellation with four satellites around L1 libration in this paper. Future work will be focused on a more complex constellation with satellites around different librations.

https://doi.org/10.1155/2017/9734164

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This study was supported by the Basic Scientific Research Fund of National Defense (no. 2016110C019), the Natural Science Foundation of Jiangsu Province (no. BK20160811), and Shanghai Deep Space Detection Technology Key Laboratory Open Funding (no. DS2016-01).

References

[1] R. H. Battin, An Introduction to the Mathematics and Methods of Astrodynamics, American Institute of Aeronautics and Astronautics, Reston, VA, USA, 1999.

[2] S. Bhaskaran, J. E. Riedel, S. P. Synnott, and T. C. Wang, "The deep space 1 autonomous navigation system: a post-flight analysis," in Proceedings of the AIAA/AAS Astrodynamics Specialist Conference, AIAA-2000-3935, Denver, CO, USA, 2000.

[3] N. Mastrodemos, D. G. Kubitschek, and S. P. Synnott, "Autonomous navigation for the deep impact mission encounter with comet tempel 1," Space Science Reviews, vol. 117, no. 1, pp. 95-121, 2005.

[4] G. S. Downs, "Interplanetary Navigation Using Pulsating Radio Sources," Tech. Rep., NASA Technical Reports, 1974, N74-34150/4.

[5] S. I. Sheikh, D. J. Pines, P. S. Ray, K. S. Wood, M. N. Lovellette, and M. T. Wolff, "Spacecraft navigation using X-ray pulsars," Journal of Guidance, Control, and Dynamics, vol. 29, no. 1, pp. 49-63, 2006.

[6] K. A. Hill, Autonomous Navigation in Libration Point Orbits [Ph.D. thesis], University of Colorado, Colorado, USA, 2007, Ph.D. dissertation.

[7] L. Zhang and B. Xu, "A universe light house--candidate architectures of the libration point satellite navigation system," Journal of Navigation, vol. 67, no. 5, pp. 737-752, 2014.

[8] L. Zhang and B. Xu, "Navigation performance of the libration point satellite navigation system in cislunar space," Journal of Navigation, vol. 68, no. 2, pp. 367-382, 2015.

[9] L. Zhang and B. Xu, "Navigation performance of the libration point satellite navigation system for future Mars exploration," Journal of Navigation, vol. 69, no. 1, pp. 41-56, 2016.

[10] Y.-J. Qian, W.-X. Jing, C.-S. Gao, and W.-S. Wei, "Autonomous orbit determination for quasi-periodic orbit about the translunar libration point," Journal of Astronautics, vol. 34, no. 5, pp. 625-633, 2013.

[11] L. Du, Z. K. Zhang, L. Yu, and S. J. Chen, "SST orbit determination of halo-lmo constellation in CRTBP," Acta Geodaetica et Cartographica Sinica, pp. 184-190, 2013.

[12] Y. T. Gao, B. Xu, and L. Zhang, "Feasibility study of autonomous orbit determination using only the crosslink range measurement for a combined navigation constellation," Chinese Journal of Aeronautics, vol. 27, no. 5, pp. 1199-1210, 2014.

[13] L. Liu and X. Y. Hou, Deep Space Probe Orbit Mechanics, Publishing House of Electronics Industry, Beijing, China, 2012.

[14] M. Olle and J. R. Pacha, "The 3D elliptic restricted three-body problem: periodic orbits which bifurcate from limiting restricted problems," Astronomy and Astrophysics, vol. 351, pp. 1149-1164, 1999.

[15] J. M. Cors, C. Pinyol, and J. Soler, "Periodic solutions in the spatial elliptic restricted three-body problem," Physica D: Nonlinear Phenomena, vol. 154, no. 3-4, pp. 195-206, 2001.

[16] J. Singh and A. Umar, "Motion in the photogravitational elliptic restricted three-body problem under an oblate primary," The Astronomical Journal, vol. 143, no. 5, article 109, 22 pages, 2012.

[17] N. Harl, K. Rajagopal, and S. N. Balakrishnan, "Neural network based modified state observer for orbit uncertainty estimation," Journal of Guidance, Control, and Dynamics, vol. 36, no. 4, pp. 1194-1209, 2013.

[18] V. Szebehely, Theory of Orbits, Academic Press, New York, USA, 1967.

[19] H. Lei, B. Xu, X. Hou, and Y. Sun, "High-order solutions of invariant manifolds associated with libration point orbits in the elliptic restricted three-body system," Celestial Mechanics and Dynamical Astronomy, vol. 117, no. 4, pp. 349-384, 2013.

10 International Journal of Aerospace Engineering [20] X. Y. Hou and L. Liu, "On motions around the collinear libration points in the elliptic restricted three-body problem," Monthly Notices of the Royal Astronomical Society, vol. 415, no. 4, pp. 3552-3560, 2011.

[21] P. Ioannou and J. Sun, Robust Adaptive Control, Prentice-Hall, Englewood Cliffs, NJ, USA, 1996.

Youtao Gao, (1) Tanran Zhao, (1) Bingyu Jin, (1) Junkang Chen, (1) and Bo Xu (2)

(1) College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing, China

(2) College of Astronomy and Space Science, Nanjing University, Nanjing, China

Correspondence should be addressed to Youtao Gao; ytgao@nuaa.edu.cn

Received 1 March 2017; Revised 11 May 2017; Accepted 24 May 2017; Published 21 June 2017

Caption: Figure 1: The barycenter synodic coordinates C-XYZ and the libration points.

Caption: Figure 2: The [L.sub.1]-centered synodic reference frame.

Caption: Figure 3: The quasiperiodic orbits around the L1 libration.

Caption: Figure 4

Caption: Figure 5

Caption: Figure 6

Caption: Figure 7
```Table 1: Initial state of the Lagrangian satellites.

Satellites           S1                    S2

x (L)        -0.009158653890369    -0.008426021835493
y (L)         0.057117380302627     0.063116578748419
z (L)        -0.005441690073741    -0.009538268525430
x (L/T)       0.071242540915778     0.077724037910168
[??] (L/T)    0.025237253422237     0.044432703977482
[??](L/T)     0.024768142966689     0.092125880718452

Satellites           S3                    S4

x (L)        -0.011852455058632    -0.012153563498421
y (L)         0.018916778218648    -0.000027358900568
z (L)         0.000050795432543    -0.001861700547480
x (L/T)       0.043377272215699     0.030872891241568
[??] (L/T)    0.002938265403217    -0.000128391119819
[??](L/T)    -0.000042323140915     0.011290787920073

Table 2: The simulation conditions.

Initial state error (m)   Observation error (m)

Scenario 1              1                        1
Scenario 2             10                        1
Scenario 3             100                       1
```