- Email: [email protected]
Attitude motion path planning on Lyapunov periodic orbits in the circular restricted three-body problem
Acta Astronautica 163 (2019) 168–180
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro
Research paper
Attitude motion path planning on Lyapunov periodic orbits in the circular restricted three-body problem
T
Mo-yao Yua, Ya-zhong Luob,∗ a b
Department of Aerospace Science and Technology, Space Engineering University, Beijing 101416, P.R. China College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, P.R. China
ARTICLE INFO
ABSTRACT
Keywords: Spacecraft attitude Path planning Three-body problem Phase plane Periodic solutions
The path of attitude motion for a rigid body in fully nonlinear Lyapunov periodic orbits is planned in the restricted three-body problem. The attitude dynamic characteristics in phase plane are analyzed, and the path constraint conditions are proposed to avoid the divergence of phase trajectories. Then, the path planning method for the three-body problem is proposed by employing the Pseudospectral method. Based on this, solving strategy of attitude-orbit coupled periodic solutions is obtained by adding event constraint conditions. Compared with previous work, the proposed method effectively avoids the divergence of attitude motion and saves a large amount of work on determining the initial guess. To verify the effectiveness of the proposed method, a large number of planning work are carried out. Besides, the paper addresses the question that how translational motions and configuration of the spacecraft affect the planned results.
1. Introduction
[15] investigated attitude motions on large nonlinear periodic orbits around Larangian points. More recently, natural periodic orbit-attitude behaviors in threebody problem have been studied. Guzzetti et al. [16,17] studied the periodic solutions of the pitch motion using Poincaré mapping, Floquet theory and differential correctors. Dayung et al. [18] explored the periodic solutions of pitch motion in the elliptic restricted three-body problem (ERTBP). However, the solutions obtained by them can only ensure the periodicity of the pitch motion, while the orientation of spacecraft may change significantly over one revolution. As denoted by Guzzetti et al. and Dayung et al., the obtainment of periodic orbit-attitude solutions is not easy. A large amount of work are carried out to find the initial guess of periodic motion because of the complex dynamics. In this paper, the critical point is to plan proper attitude maneuver path that can keep moving without divergence. An expected orientation is artificially selected in advance and the planned path is expected to move around it as close as possible. The planning method is proposed based on the Pseudospectral method and the dynamic analysis on phase plane. Then, boundary conditions, convergent conditions and path constraints are deduced. In contrast to previous studies [13–16] active expected factors are considered instead of passive natural motion. At the same time, divergence of pitch angle can be significantly reduced. Base on the proposed method, coupled attitude-orbit periodic solutions
In the past three decades, the completion of several libration orbit missions, such as the International Sun-Earth Explorer-3 [1], Genesis mission [2] and some others [3], has proved that the Lagrangian points are of great use, which attracts much attention in the exploration of three-body problem. However, because of the sensitive environment resulting from the coupled-gravitational field, spacecraft exhibit complicated behaviors. Thus, understanding, predicting and finally planning the attitude motion in the complex dynamical environment are increasingly significant. The circular restricted three-body problem (CR3BP) has been studied extensively in recent years [4,5]. However, the attitude dynamic problem in this regime has not been fully explored. In contrast to the two-body problem, spacecraft at a Lagrangian point or moving around it show complicated dynamic characteristics. The gravity gradient torques by the two primary bodies differ as the attitude or position of the spacecraft changes. Thus, the attitude dynamics of spacecraft at the Lagrangian points and their vicinity show new features compared to the two-body system. The earliest investigations by Kane et al. [6,7], Robinson [8,9], Abad [10] and Brucker et al. [11] focused on spacecraft that are artificially fixed at the Lagrangian points. Then, the effects of gravity torques along small linear periodic orbits are examined by Wong et al. [12], while Lara [13], Yunhe et al. [14] and Knutson et al.
∗
Corresponding author. E-mail addresses: [email protected] (M.-y. Yu), [email protected] (Y.-z. Luo).
https://doi.org/10.1016/j.actaastro.2019.08.013 Received 29 July 2018; Received in revised form 11 August 2019; Accepted 13 August 2019 Available online 21 August 2019 0094-5765/ © 2019 IAA. Published by Elsevier Ltd. All rights reserved.
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
where 1
=
(x + µ ) y z , µ1 = , 1 = r1 r1 r1 (x + µ 1) y , µ2 = , 2 = r2 r2
2
=
z r2
The orbital motion of the spacecraft under the gravity by P1 and P2 is expressed in the primary frame as
y¨ = Fig. 1. Geometry, notation, and frames of the problem.
z¨ =
µ
1
x¨ = 2y + x
r13
2x + y 1
µ r13
z
(x + µ ) 1
µ
r13 µ z r23
y
µ (x + µ r23 µ y r23
1)
Since the late 1960s, several types of periodic motions around Lagrangian points have been studied under the background of CR3BP, such as the Lyapunov orbits, Halo orbits and Lissajous orbits. These periodic orbits are usually computed through analytic method or numerical method. Early researches mostly use theoretical analytic and analytic calculation [19,20]. With the development of digital computers, the numerical solutions of periodic orbits attract more and more attention [21–23]. Compared with analytic method, the obtainment of numerical solutions is easier. In this study, the out-plane motion of the spacecraft is ignored, only the periodic motion in the plane of the two primary bodies is considered, employing the numerical method. The planar Lyapunov orbit are obtained by solving the first two equations of Eq. (3) as a Two-Point Boundary Value Problem (TPBVP) [24]. Considering the symmetry of the periodic orbit about axis l1, the problem can be transformed into an Initial Value Problem (IVP) [25]. The initial position is usually chosen on axis l1, thus the initial parameters for Lyapunov orbit are
around the expected orientation can be obtained by adding event constraint conditions. Large amount of efforts can be saved on finding the initial value. Finally, a large number of planning work are carried out to verify the effectiveness of the proposed method. 2. Model The finite-sized rigid body spacecraft is considered in the present paper. In Fig. 1, P1 is the more massive attractor with mass m1, and P2 is the second attractor with mass m2. The mass ratio between P1 and P2 is µ = m2 /(m1 + m2 ) . The distance between P1 and P2 is denoted as d. The dynamics of the spacecraft are modeled under the assumption of planar circular restricted three-body problem (CR3BP), such that: 1) the size and the mass of the spacecraft can be ignored compared with the twoprimary bodies; and 2) P1 and P2 move on a circular orbit at an angular velocity Ω around O, which is the center of mass of P1 and P2. Before specifically developing the work, it is necessary to supply a general framework. Let a primary frame l, rotating with the primaries P1 and P2, be fixed at O. The three unit component vectors are denoted as l1, l2, and l3. l1 is aligned from P1 to P2, l3 is normal to the orbital plane of the primary bodies, and l2 completes a right-handed set of orthogonal axes with l1 and l3. Define the local frame b of the spacecraft to be fixed in the rigid body and attached to the center of mass o. The three unit component vectors denoted as b1, b2, and b3 are along the principal axes of inertia of the spacecraft. The two frames are related by yaw angle , pitch angle , and roll angle , with a rotation sequence of 3-2-1 from the primary frame to the local frame. For a finite-sized rigid spacecraft, a complete description of motion includes the orbital and attitude dynamics. Because of the significantly different scales of variables involved in these two kinds of motion, computational issues may arise. To overcome this problem, the units of mass, length and time are chosen as the sum of mass, the distance and the angular velocity of the two-primary system, respectively, such that m1+m2=1, d=1, and Ω =1. Obviously, the gravitational constant then has the value of unity.
[x y x y ]T = [x 0 0 0 y0 ]T
2.2. The attitude motion Based on the assumption that the spacecraft is influenced only by P1
Denote the position of the spacecraft as [x, y, z] in the primary frame, which is a position vector with respect to the barycenter. Then vectors pointing to the spacecraft from P1 and P2 can be expressed as
r1 = (x + µ ) l1 + yl2 + zl3 1) l1 + yl2 + zl3
(1)
Parameter µ in Eq. (1) is the distance from P1 to O, which is actually equal to the mass ratio after non-dimension. The unit coefficients of vectors r1 and r2 can be expressed as
r1 = r2 =
1 l1 2 l1
+ µ1 l2 + + µ 2 l2 +
1 l3 2 l3
(4)
For a given x0, the orbit is obtained by integrating Eq. (3) until the final y and integral time are changeless within the given tolerance. Fig. 2 shows two L1 planar Lyapunov orbits in the Earth-Moon system. The orbits in Fig. 2 span an amplitude Ay (maximum coordinate value in the l2 direction) of 0.2862 and 0.3331 after non-dimension, respectively.
2.1. The orbital motion
r2 = (x + µ
(3)
Fig. 2. Two planar Lyapunov orbits.
(2) 169
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 3. Phase trajectories of L1 Lagrangian point in the Earth-Moon three-body system.
and P2, while other environmental torques are treated as perturbations. For a rigid spacecraft subjected to gravity gradient torques, the equation of attitude motion of in the CR3BP is described by Euler's equation as
I
+
×I
=
single first-order equation by substituting Eq. (1) and Eq. (6) into Eq. (5)
1¨ =
(5)
g
Here, I is the spacecraft's inertia matrix in diagonal form, denoted as I = diag [I1, I2, I3]. denotes the body frame's angular velocity with respect to the inertial space. The gravity gradient torque g due to P1 and P2 is given as g
=
3m1 3m r1 × Ir1 + 5 2 r2 × Ir2 r15 r2
A=
sin cos 0
0 0 1
2 i=1
3mi ( ri3
2 2 i
Here the coefficient
=
I2
µi2 )sin 2
+ i=1
3m i ri3
i µi
cos 2
(10)
represents the inertia ratio
I1 I3
(11)
The subsequent analysis is carried out based on the attitude equation (10) and the orbital equation (3).
(6)
Eq. (3) and Eq. (5) are solved together to obtain the orbit-attitude solutions of the spacecraft. In this study, rotations and displacements outside the plane P1 and P2 are ignored, that means, z = 0 , 1 = 2 = 0 , and the rotation matrix A is a function of the pitch angle
cos sin 0
1 2
3. Methodology In this section, the attitude dynamic characteristics in phase plane are firstly analyzed. Based on it, boundary conditions, convergent conditions and path constraints are deduced. In the final subsection, the Pseudospectral method for path planning is described.
(7)
Additionally,
l1 = [ cos l2 = [ sin
sin cos
l3 = [0 0 1]T
3.1. Dynamic analysis on phase plane
0]T 0]T
To get the phase trajectories, integrate Eq. (10) over time. Then the following equation can be obtained
(8)
The angular velocity in the body frame with respect to the inertial frame is
=
+
2
+
here,
(9)
Consequently, the description of attitude behavior reduces to a 170
sin(2 + ) = E
(12)
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 4. The Phase Angle
=
and d0 in the Earth-Moon Space.
curve in Fig. 3(b) is greater than that in Fig. 3(a), and it is the same for Fig. 3(c) and (d). According to Eq. (12) and Eq. (13), the amplitude of the critical curve can be expressed as
c12 + c22 2
1 2
c1 =
i=1 2
c2 = i=1
sin =
c1
3mi ( ri3
3mi ri3 , cos
2 i
µi2 )
d0 =
i µi
=
c2
(13)
=
E
(
)Stablemotion around the equilibrium
Critical stable(sinecurvespassingthrough , )divergent motion(wavycurves)
4
(17)
2
0. 2192
(18)
These two solutions correspond to the center of the two special regions in Fig. 4, respectively. As the amplitude d0 approaches zero, the stable regions reduce to a narrow ellipse and finally turn out to be an equilibrium point. Then the closed curves vanish, without any lines passing through the axis = 0 . Thus, all the attitude state including the equilibrium points become unstable, and a diverging pitch motion occurs. Actually, substitute = 0.2 and = 0 into Eq. (12), the following figures can be drawn. Fig. 5(a) shows the phase trajectories near the center of the special region. The small amplitude of the critical curve in Fig. 5(a) limits the stable region to a narrow one. Fig. 5(b) shows the central phase trajectories of the special region. It can be noted that the stable region vanished. Actually, the red lines in Fig. 5(b) correspond to infinite discrete points on axis = 0 . The dynamic characteristics of these points are all balanced. Different from the equilibrium points in Fig. 3, the phase trajectories around these points turn into straight lines instead of closed curves, sine curves or wavy curves. That means the
point(closed curve) E0 = E (
or
x2 = 0. 9373, y2 =
Thestableequlibriumpoint ,
2
x1 = 0. 9373, y1 = 0. 2192
To sum up, the value range of E and the corresponding characteristics on the phase plane are as follows.
Emin =
3 4
d0 determines the size of stable region. The bigger d0 is, the larger the stable region is. Fig. 4 shows the distribution of and d0 in the EarthMoon space. Two regions emerge on the upside and downside of the moon in both Fig. 4(a) and (b). The phase angle changes sharply from −90° to 90° and the amplitude of the critical curve d0 reaches zero in these two regions. Actually, when solving the equation = 0 , two solutions about x and y can be obtained for the Earth-Moon system.
(14)
sin(2 + )) =
(16)
2
If the position of the spacecraft is fixed, E0 is simply decided by . Thus, a larger value of produces a larger amplitude of sine curve. Generally speaking, for a fixed position in the three-body space, four equilibrium points of the pitch angle motion can be found. Their characteristics of stability is decided by the sign of the inertia ratio , and the size of the stable regions is determined by the value of , which ranges from −1 to +1. The above research obtains the characteristics of phase trajectories at a fix position. When it comes to the whole three-body space, it is important to have a general understand of the distribution of equilibrium angle and the size of stable region. According to Eq. (12) and Eq. (16), the phase angle determines the position of equilibrium angle. The equilibrium angle can be expressed as
Note that, E is determined by the position and attitude of the spacecraft, ranging from to positive infinity. For a fix position in the three-body plane, the value of is constant. Then families of phase path can be drawn by changing the value of E. Fig. 3 shows the phase trajectories at L1 Lagrangian point for different inertia ratios. Four equilibrium points can be noted on the phase plane in Fig. 3, with a difference of 90° between neighbors, denoted by red points. The shadow regions represent stable attitude motions. When the attitude state of the spacecraft falls into these regions, it will move on along the closed curve. Thus, the attitude angle will vibrate periodically around the center equilibrium point. However, when the attitude state is beyond these regions, it will change along the wavy curves. In this case, the attitude angle will diverge. If the attitude state is above the shadow region, the angle velocity is positive, and the attitude angle will grow to positive infinity and vice versa. Obviously, the two saddle points are unstable while the other two points are stable. The stable points correspond to the minimum value of E, while the saddle points correspond to a critical value E0. With the increase of E, divergence occurs on the pitch angle motion after it reaches E0. The phase trajectories no longer form a closed curve, turning into a wavy line. Actually, when E is greater than E0, cannot reach the zero point, which means that E is sin(2 + ) . Obviously, the crilarger than the maximum value of tical value of E0 is
E0 = max(
2E0 =
axis) (15)
Compared with Fig. 3(a) and (c), Fig. 3(b) and(d), it can be discovered that the stable equilibrium points between them take a difference of π/2, which can be explained by: a) If > 0 , E0 = and = /4 /2 . b) If < 0 , E0 = /2 . According to and = /4 Eq. (13), for a fixed position, is constant. Thus, changing the sign of leads to a movement of π/2 for stable equilibrium points. Then considering Fig. 3(a) and (b), the amplitude of the critical 171
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 5. Phase trajectories in the two special regions.
Fig. 6. Phase trajectory for Ay = 0.2862 ( = 0.6 ).
balance is easily to be disturbed, and becomes divergent quickly. Actually, the balanced points at the center of the special region are much more unstable than the saddle points discussed before. However, almost
all the common periodic orbits, including Lyapunov orbits and Halo orbits around L1 and L2 Lagrangion points, transit through these regions. According to the results of Knutson et al. [15], spacecraft can 172
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 7. Phase trajectory for Ay = 0.3331 ( = 0.6 ).
hardly remain close to its initial orientation on planar Lyapunov orbit except for certain spacecraft topology and orbit amplitude. Simulations were carried out by Knutson et al. [15], and it was concluded that the spacecraft can remain close to its original attitude ( = = 0 ) on orbits Ay = 0.2862. However, when it comes to orbit Ay = 0.3331, a great change takes place in pitch motion, nearly 180°. Figs. 6 and 7 shows the attitude motions on these two orbits over one revolution in the phase plane. Different figures correspond to different moment in one revolution, as shown in the subtitle. The red curves are the phase trajectories of the spacecraft from the beginning to the
current time. The blue points correspond to the current attitude state, that means the horizontal and vertical coordinates of this point represent the current attitude angle and the attitude velocity, respectively. The black curves are the phase trajectories of the current orbital position. Take Fig. 6 (a) as an example, the horizontal and vertical coordinates of the blue point are the attitude state at 0.34 revolution and the phase trajectory corresponds to the orbital position at 0.34 revolution. It is shown in Fig. 6 that the phase trajectory change a lot from 0.4 revolution to 0.6 revolution, which actually corresponds to the two 173
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 8. The planning phase trajectories for Ay = 0.3331 ( = 0.6 ).
through half a revolution, which directly leads to the divergence of the pitch motion. Therefor, when the phase trajectory moves beyond the stable region, divergence is most likely to take place. 3.2. The path planning schedule Equations of Motion. In this study, we are focused mainly on the global orbit-attitude dynamics of the pitch motion in the circular restricted three-body problem. Thus, the attitude equation (10) and the orbital equation (3) derived above are solved together to achieve the orbit-attitude solutions of the spacecraft. Objective Function. It is usually essential for spacecraft to maintain a certain orientation. Thus, an expected pitch angle is selected, denoted as ˜ (t ) . Then the objective is to plan a maneuver path near the expected orientation as close as possible. In other words, the maximum deviation from the expected orientation should be the smallest. Thus, the objective function can be expressed as
Fig. 9. Spacecraft Configuration example.
min J = max(
(t )
~
(19)
(t ) )
Constraints. The constraints contain boundary conditions and path constraints. In some special cases, the event constraints are also employed. In this study, no artificial control is introduced. The attitude motion is decided simply by the initial attitude and position. The path planning work is actually to find the proper initial state and to achieve the greatest satisfaction. So no boundary conditions of the attitude motion are employed here. The boundary conditions of the orbital motion can be expressed as (20)
x (t 0) = x 0 , y (t 0) = 0, x (t 0 ) = 0, y (t0) = y0
To ensure the stability, the path constraints are set to keep the phase trajectory inside the stable region of the local position. Thus, the path constraints are denoted as
Fig. 10. The Pitch Angle Response on L1 Lyapunov family over one revolution. (Earth-Moon system).
E (t ) =
(t ) 2 +
(t )sin(2 (t ) + (t )) < E0 (t ) =
(t )
(21)
Additionally, periodic solutions can be obtained by adding event constraints
special regions in Fig. 4. During this period, the phase angle changes by about 180°and the amplitude of the critical curve varies from smaller than 1 to larger than 10. The same characteristics can be found in Fig. 7. However, the essential difference between Figs. 6 and 7 is the stability of the attitude motion. The motion on orbit Ay = 0.2862 is always within the scope of stable region, as shown by the black curves. On the contrary, the motion in Fig. 7 enters the unstable region after passing
(t 0) =
(t f )
(t 0) =
(t f )
(22)
Here the subscript 0 means the initial time and f represents the final time. 174
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 11. The Pitch Angle Response on L1 Lyapunov family over multi revolutions. (Earth-Moon system).
3.3. Pseudospectral method
the algebraic equations at LGL nodes and the optimal control problem is transformed into a nonlinear programming (NLP). Then, the NLP can be solved by sequential quadratic programming algorithm. In this paper, the GPOPS is used, which is a PSM Matlab tool box design by Rao [26].
The Legendre pseudospectral method is employed here. The states are estimated with Legendre polynomials Pi that is orthogonal on the interval [-1, 1]. The nth degree Legendre polynomial is
1
Pn (t ) =
2nn!
dn 2 (t dt n
1)n
4. Results
(23)
4.1. Method verification
The states are interpolated at pre-selected node points, and the values at these points become the optimal variables. The dynamic equation (9) and (3) and the path constraints are enforced at the Legendre-Gauss-Labatto (LGL) points. The constraints for the entire interval are assumed to be satisfied when they are satisfied as the node points. Define x (t ) = [x (t ) y (t ) x (t ) y (t ) (t ) (t )]T to be the vector of all states for this problem and let x (t ) = f (x (t ), v~ (t )) represent the corresponding right-hand sides of the dynamic Eq. (9) and (3). For any given time interval [t 0, t f ], the transformation t : [ t0, tf ] = [ 1,1] [t 0, t f ] is given by
(t f
t (t ) =
t 0 ) t + (t f + t 0 )
To begin the path planning work in the three-body problem, an example is employed to verify the effectiveness of the proposed method. The divergence motion on orbit Ay = 0.3331 is planned, which is shown in Fig. 7. To make a comparison, the expected pitch motion is chosen as ˜ (t ) = 0 . Fig. 8 shows the planned results. The planned initial attitude angle is (t 0 ) = 7.15 , just with a little deviation from (t 0 ) = 0 which can be ignored. After planning, the phase trajectory locates in the stable regions over the whole revolution. The maximum deviation from the expected pitch angle is only about 13°. Divergence does not occur here, which demonstrates that the proposed planning method is effective for this case. Additionally, the planned phase trajectories seem to be symmetry, which may indicate a better dynamic characteristic.
(24)
2
Eq. (23) can be used to shift the LGL points to the desired interval. Using the LGL points, the state variables are estimated by
4.2. Mapping the pitch angle response
N
x N (t ( t )) =
x (t )
xi
i (t
)
(25)
i=0
where xi are coefficients, satisfies i (t
i (t
) is the Nth Lagrange polynomials that
1, if i = k 0, if i k
)=
This subsection is focused on planning the attitude motion path along different Lyapunov periodic orbits. As the amplitude of periodic orbits increases, the fully nonlinear dynamics over one revolution become more complicated, and the phase angle and d0 also change greatly. As discussed above, Knutson et al. [15] have carried out a large amount of simulations to investigate the natural attitude response along these orbits. The initial condition was set to be (t ) = (t ) = 0 by them. In the present work, the L1 family is selected again, two orbits of which are displayed in Fig. 2. The planned results will be shown in a map to compare with the natural response by Knutson et al. The variables in the map are chosen as the amplitude of each periodic orbit and the configuration of the spacecraft. The configuration of the spacecraft is described by the inertia ratio , as defined by Eq. (11). Take the model in Fig. 9 as an example to illustrate the relation between the configurations and the inertia ratio . The spacecraft is regarded as an approximate cuboid, the three body vectors are along the principal axes of inertia of the spacecraft and the three side lengths are denoted as a, b and c, respectively. Then Eq. (29) can be obtained. It is observed that when ignoring the rotations and displacements outside of the orbital plane of P1 and P2, the spread of the spacecraft in this direction makes no effects on the attitude motion. Instead of setting initial states, the expected attitude motion is
(26)
The derivatives of the state variables at the LGL nodes can be computed by
x N (t ( t )) =
x (t ( t ))
N
2 (t f
t0 )
xi i=0
i (t
)=
N
2 (t f
t0 )
Dki x i i=0
(27)
The (N + 1) × (N + 1) differential matrix D associated with the Legendre polynomials is PN (tk ) PN (tk )(t f ti )
Dki =
N (N + 1) 4 N (N + 1) 4
0
if k
i
if k = i=0 if k = i=N otherwise
(28)
Thus, the differential equation of the control model is estimated by 175
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 12. Multi-Period Simulation results. (Ay = 0.2037, κ = 0.1).
selected to be ˜ (t ) = 0 for the sake of comparison. Maintaining this expected attitude, different simulations are completed by varying the amplitude of the orbit and the inertia ratio . Fig. 10 is a map of the global deviation from the expected attitude over one revolution. Every point in the map corresponds to a planning work over one revolution, then the maximum bias of the pitch angle from the expected motion is recorded. Different colors in Fig. 10 represent the value of the bias. The computation time for one planning task is averagely 20s with one processor (8 GB RAM).
If = 1, a=0 If 1 < < 0, a < b If = 0, a=b a>b If 0 < < 1, b=0 If = 1,
(29)
On the whole, the planned pitch angle over one revolution on the L1 Lyapunov family maintains closer to the expected orientation than the natural motion with the expected attitude being its initial condition. 176
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 13. Multi-Period Simulation results. (Ay = 0.4556, κ = −0.6).
Besides, all the phase trajectories of the planned results fall into the stable regions on the phase plane. Although the pitch angle grows larger with the increase of Ay, bias larger than 70 deg has not been found within the scope of this study. Compared with the natural motion in Ref. [13], the proposed planning method works well for the fully nonlinear attitude-orbit coupled system. Additionally, a dark horizontal band emerges from the map for a narrow set of Lyapunov orbits, which is consistent with the feature discovered by Knutson. Fig. 11 gives the maps for 2 revolutions and 4 revolutions, respectively. As is mentioned above, the computation time for a one-
revolution planning work costs about 20s. When it comes to multiperiod planning work, a tremendous growth in planning time occurs since the LGL points and the Legendre polynomials increase with the periodicity. The overall characteristics of distribution in Fig. 11 are similar to that in Fig. 10. When the spacecraft is axisymmetric about b3, it can always maintain = = 0 . With the growth of amplitude of the orbit, the pitch angle diverges more greatly from the expected orientation, especially for 4-revolution planning results. The largest angle for 4-revolution reaches 90° in Fig. 11(b). It is worth noting that the dark horizontal band turns less obvious as the planning period 177
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 14. The phase trajectories of periodic solutions on orbit Ay = 0.1817. ( ˜ (t ) = 0 ).
Fig. 15. The phase trajectories of periodic solutions on orbit Ay = 0.3044. ( ˜ (t ) = 0 ).
Fig. 16. The phase trajectories of periodic solutions on orbit Ay = 0.1817. ( ˜ (t ) = /2 ).
Fig. 17. The phase trajectories of periodic solutions on orbit Ay = 0.3044. ( ˜ (t ) = /2 ).
Fig. 18. The phase trajectories of periodic solutions on orbit Ay = 0.1817. ( ˜ (t ) =
178
/2 ).
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 19. The phase trajectories of periodic solutions on orbit Ay = 0.3044. ( ˜ (t ) =
increases. It can be predicted that this band will disappear with the increase of periodicity. Since the multi-period planning work normally takes much time, multi-period simulations are carried out with initial conditions determined by one-revolution planned results, aiming to see whether the one-period planned results work well for multi periods. Two couples of orbits and spacecraft configurations are selected, one of which falls into the regions less than 10° in Fig. 10, and the other one is in the regions more than 50°. The simulation results are shown in Figs. 12 and 13. The multi-period simulations of the selected orbits show good characteristics on maintaining the attitude motion around the expected orientation. In Fig. 12, the phase trajectories keep a circumferential motion in an adjacent area. The pitch angle shows a waving curve less than 10°. Actually, this stable motion can keep going for a much longer simulation time. In Fig. 13, the phase trajectories show a similar-saddle shape. Within the time frame of the study (4 revolution), no divergence occurs. Although the selected orbits cannot represent all the situations, the superiority of the proposed planning method is outstanding when compared with the none-planning simulations.
/2 ).
feature in these phase trajectories is that they are all closed curves, which guarantees the periodicity. Compared with other methods that aim to find periodic solutions, a large amount of efforts on determining the initial guess can be saved by bringing in the expected orientation. 5. Conclusion The path of attitude motion for a rigid body in fully nonlinear Lyapunov periodic orbits is planned in the background of restricted three-body problem. The attitude dynamic characteristics in phase plane are deeply analyzed, and the path constraint conditions are proposed to avoid the divergence of phase trajectories. Then, a path planning method in three-body problem is proposed by employing the Pseudospectral method. Based on this, solving strategy of attitude-orbit coupled periodic solutions is obtained by adding event constraint conditions. A large number of simulations are carried out to explore the attitude behavior across the planar Lyapunov families by changing the configuration of the spacecraft. Results are summarized in maps, which address the question that how the translational motions and the configuration of the satellite affect the planned results. Finally, multiperiod results and attitude-orbit periodic solutions are conducted to make comparison with previous researches. It is shown that the proposed method effectively avoids the divergence of attitude motion, and the bias of pitch angle can be significantly reduced. Besides, a large amount of efforts on determining the initial guess can be saved.
4.3. Periodic solutions As is mentioned before, periodic attitude-orbit solutions can be obtained by adding event constraints described in Eq. (22). It is more complicated that for a selected reference orbit, there may be several periodic solutions, corresponding to different initial pitch angles, which have been discovered by Dayung and Rodney [18]. In their study, not only one-period attitude-orbit solutions but also several multi-period solutions come out for one orbit. In the present study, the expected orientation is also brought in, aiming to select the most satisfied periodic solutions. Based on the existing planning method, event constraints are introduced to guarantee that the pitch angle and angular velocity at the initial time are equal to that at the final time in a revolution. Thus, periodic solutions can be obtained by the planning work. It is worth noting that the attitude-orbit periodic solutions raise a high requirement on the accuracy of periodic orbit. The pre-set tolerance of the iteration to get Lyapunov periodic orbits is 10−5. Two orbits are chosen with an amplitude Ay of 0.1817 and 0.3044 to show the planning periodic solutions in this paper. Three cases are planned with different expected orientations. Case 1: ˜ (t ) = 0 Case 2: ˜ (t ) = /2 /2 Case 3: ˜ (t ) = Figs. 14–19 show the phase trajectories of the planned periodic solutions. Three expected orientations are chosen as ˜ (t ) = 0 , ˜ (t ) = /2 and ˜ (t ) = /2 , corresponding to Case 1, Case 2 and Case 3, respectively. It can be found that the phase trajectories show great difference in shape and size for different orbits, inertia ratios and expected orientations. Generally speaking, the planned trajectories of pitch angle are all moving around the given expected orientation. Certainly, the periodic solutions on one orbit is limited. Not every given expected orientation has a unique corresponding curve. The common
6. Further investigation The planned method in this paper aims to explore the characteristic of attitude motion under the background of CR3BP. But when a certain specific mission is considered, the solutions may be unaccepted. For example, the 70 deg bias over one revolution is not allowed for an observation mission. Thus, the planned method should be improved to consider control [28]. Besides, applying the planning method in the Halo orbit with three dimensions [29] is worthy of further exploration. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 11572345). References [1] R.W. Farquhar, The flight of ISEE-3/ICE: origins, mission history, and a legacy, J. Astronaut. Sci. 49 (1) (2001) 23–74. [2] M.W. Lo, B.G. Williams, W.E. Bollman, D. Han, Y. Hahn, J.L. Bell, B. Barden, Genesis mission design, J. Astronaut. Sci. 49 (1) (2001) 169–184. [3] D.W. Dunham, R.W. Farquhar, G. Gómez, M.W. Lo, J.J. Masdemont (Eds.), “Libration Point Mission, 1978–2002,” Libration Point Orbits and Applications: Proceedings of the Conference, World Scientific Publ. Company, Aiguablava, June 2003. [4] L. Bucci, M. Lavagna, D. Guzzetti, Periodic orbit-attitude solutions along planar orbits in a perturbed circular restricted three-body problem for the Earth-Moon System, Acta Astronaut. 99 (2018) 152–162. [5] D. Guzzetti, K.C. Howell, Attitude dynamics in the circular restricted three-body problem, Astrodynamics 2 (No. 2) (2018) 87–119.
179
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo [6] T.R. Kane, E.L. Marsh, Attitude stability of a symmetric satellite at the equilibrium points in the restricted three-body problem, Celest. Mech. Dyn. Astron. 4 (No. 1) (1971) 78–90. [7] T.R. Kane, P.M. Barba, Attitude stability of a spinning satellite in an elliptic orbit, J. Appl. Mech. 33 (No. 2) (1966) 402–405. [8] W.J. Robinson, Attitude stability of a rigid body placed at an equilibrium point in the restricted problem of three bodies, Celest. Mech. Dyn. Astron. 10 (No. 1) (1974) 17–33. [9] W.J. Robinson, The Restricted problem of three bodies with rigid dumb-bell satellite, Celest. Mech. Dyn. Astron. 8 (2) (1973) 323–330. [10] A. Abad, M. Arribas, A. Elipe, On the attitude of a spacecraft near a Lagrangian point, Astron. Inst. Czechoslovakia Bull. 40 (5) (Sept. 1989) 302–307. [11] E. Brucker, P. Gurfil, Analysis of gravity-gradient perturbed rotational dynamics at the collinear Lagrange points, J. Astronaut. Sci. 55 (3) (2007) 271–291. [12] B. Wong, R. Patil, A. Misra, Attitude dynamics of rigid bodies in the vicinity of Lagrangian points, J. Guid. Control Dyn. 31 (1) (2008) 252–256. [13] M. Lara, J. Peláez, C. Bombardelli, F.R. Lucas, M. Sanjurjo-Rivo, D. Curreli, E.C. Lorenzini, D.J. Scheeres, Dynamic stabilization of L2 periodic orbits using attitude-orbit coupling effects, J. Aerosp. Eng. 4 (1) (2012) 73–82. [14] Y.H. Meng, R. Hao, Q.F. Chen, Attitude stability analysis of a dual-spin spacecraft in halo orbits, Acta Astronaut. 99 (2014) 318–329. [15] A.J. Knutson, D. Guzzetti, K.C. Howell, Attitude responses in coupled orbit-attitude dynamical model in Earth–moon Lyapunov orbits, J. Guid. Control Dyn. 38 (7) (2015) 1264–1273. [16] D. Guzzetti, K.C. Howell, “Coupled Orbit-Attitude Dynamics in the Three-Body Problem: A Family of Orbit-Attitude Periodic Solutions,” AIAA/AAS Astrodynamics Specialist Conference, Aug. 2014 AIAA Paper 2014-4100.
[17] D. Guzzetti, K.C. Howell, Natural periodic orbit-attitude behaviors for rigid bodies in three-body periodic orbits, Acta Astronaut. 130 (2017) 97–113. [18] D. Koh, R.L. Anderson, Periodic orbit-attitude solutions in the planar elliptic restricted three-body problem, J. Guid. Control Dyn. 41 (3) (2018) 644–657. [19] H. Poincaré, New Methods of Celestial Mechanics. History of Modern Physics and Astronomy, American Institute of Physics, 1892. [20] D.L. Richardson, Analytical construction of periodic orbits about the collinear points, Celestial Mech. 22 (3) (1980) 241–253. [21] J. Breakwell, J. Brown, The “Halo” family of 3-dimensional periodic orbits in the Earth-Moon restricted 3-body problem, Celestial Mech. 20 (1979) 389–404. [22] M. Xu, T. Tan, S. Xu, Research on the transfers to Halo orbits from the view of invariant manifolds, Sci. China Phys. Mech. Astron. 55 (4) (2012) 671–683. [23] T. Luo, M. Xu, Y.F. Dong, Natural formation flying on quasi-halo orbits in the photogravitational circular restricted three-body problem, Acta Astronaut. 149 (2018) 35–46. [24] P.B. Bailey, L.F. Shampine, P.E. Waltman, Nonlinear Two Point Boundary Value Problems, Academic press, 1968. [25] J.D. Lambert, Numerical Methods for Ordinary Differential Systems: the Initial Value Problem, John Wiley & Sons, Inc., 1991. [26] A.V. Rao, D.A. Benson, Algortihm 902: GPOPS, a Matlab software for solving multiple-phase optiaml control using the gauss pseudospectral method, ACM Trans. Math Software 37 (2) (2010) 1–39. [28] M. Xu, Y.H. Jia, S.J. Xu, Attitude description and dynamics modeling for Halo-orbit spacecraft, J. Beijing Univ. Aeronaut. Astronaut. 33 (10) (2007) 1166–1169. [29] M. Xu, Y.Y. Liang, T. Luo, X. Sun, Attitude pointing schemes and spacecraft configurations for libration-point-orbit spacecraft, Aero. Sci. Technol. 65 (2017) 1–8.
180
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro
Research paper
Attitude motion path planning on Lyapunov periodic orbits in the circular restricted three-body problem
T
Mo-yao Yua, Ya-zhong Luob,∗ a b
Department of Aerospace Science and Technology, Space Engineering University, Beijing 101416, P.R. China College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, P.R. China
ARTICLE INFO
ABSTRACT
Keywords: Spacecraft attitude Path planning Three-body problem Phase plane Periodic solutions
The path of attitude motion for a rigid body in fully nonlinear Lyapunov periodic orbits is planned in the restricted three-body problem. The attitude dynamic characteristics in phase plane are analyzed, and the path constraint conditions are proposed to avoid the divergence of phase trajectories. Then, the path planning method for the three-body problem is proposed by employing the Pseudospectral method. Based on this, solving strategy of attitude-orbit coupled periodic solutions is obtained by adding event constraint conditions. Compared with previous work, the proposed method effectively avoids the divergence of attitude motion and saves a large amount of work on determining the initial guess. To verify the effectiveness of the proposed method, a large number of planning work are carried out. Besides, the paper addresses the question that how translational motions and configuration of the spacecraft affect the planned results.
1. Introduction
[15] investigated attitude motions on large nonlinear periodic orbits around Larangian points. More recently, natural periodic orbit-attitude behaviors in threebody problem have been studied. Guzzetti et al. [16,17] studied the periodic solutions of the pitch motion using Poincaré mapping, Floquet theory and differential correctors. Dayung et al. [18] explored the periodic solutions of pitch motion in the elliptic restricted three-body problem (ERTBP). However, the solutions obtained by them can only ensure the periodicity of the pitch motion, while the orientation of spacecraft may change significantly over one revolution. As denoted by Guzzetti et al. and Dayung et al., the obtainment of periodic orbit-attitude solutions is not easy. A large amount of work are carried out to find the initial guess of periodic motion because of the complex dynamics. In this paper, the critical point is to plan proper attitude maneuver path that can keep moving without divergence. An expected orientation is artificially selected in advance and the planned path is expected to move around it as close as possible. The planning method is proposed based on the Pseudospectral method and the dynamic analysis on phase plane. Then, boundary conditions, convergent conditions and path constraints are deduced. In contrast to previous studies [13–16] active expected factors are considered instead of passive natural motion. At the same time, divergence of pitch angle can be significantly reduced. Base on the proposed method, coupled attitude-orbit periodic solutions
In the past three decades, the completion of several libration orbit missions, such as the International Sun-Earth Explorer-3 [1], Genesis mission [2] and some others [3], has proved that the Lagrangian points are of great use, which attracts much attention in the exploration of three-body problem. However, because of the sensitive environment resulting from the coupled-gravitational field, spacecraft exhibit complicated behaviors. Thus, understanding, predicting and finally planning the attitude motion in the complex dynamical environment are increasingly significant. The circular restricted three-body problem (CR3BP) has been studied extensively in recent years [4,5]. However, the attitude dynamic problem in this regime has not been fully explored. In contrast to the two-body problem, spacecraft at a Lagrangian point or moving around it show complicated dynamic characteristics. The gravity gradient torques by the two primary bodies differ as the attitude or position of the spacecraft changes. Thus, the attitude dynamics of spacecraft at the Lagrangian points and their vicinity show new features compared to the two-body system. The earliest investigations by Kane et al. [6,7], Robinson [8,9], Abad [10] and Brucker et al. [11] focused on spacecraft that are artificially fixed at the Lagrangian points. Then, the effects of gravity torques along small linear periodic orbits are examined by Wong et al. [12], while Lara [13], Yunhe et al. [14] and Knutson et al.
∗
Corresponding author. E-mail addresses: [email protected] (M.-y. Yu), [email protected] (Y.-z. Luo).
https://doi.org/10.1016/j.actaastro.2019.08.013 Received 29 July 2018; Received in revised form 11 August 2019; Accepted 13 August 2019 Available online 21 August 2019 0094-5765/ © 2019 IAA. Published by Elsevier Ltd. All rights reserved.
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
where 1
=
(x + µ ) y z , µ1 = , 1 = r1 r1 r1 (x + µ 1) y , µ2 = , 2 = r2 r2
2
=
z r2
The orbital motion of the spacecraft under the gravity by P1 and P2 is expressed in the primary frame as
y¨ = Fig. 1. Geometry, notation, and frames of the problem.
z¨ =
µ
1
x¨ = 2y + x
r13
2x + y 1
µ r13
z
(x + µ ) 1
µ
r13 µ z r23
y
µ (x + µ r23 µ y r23
1)
Since the late 1960s, several types of periodic motions around Lagrangian points have been studied under the background of CR3BP, such as the Lyapunov orbits, Halo orbits and Lissajous orbits. These periodic orbits are usually computed through analytic method or numerical method. Early researches mostly use theoretical analytic and analytic calculation [19,20]. With the development of digital computers, the numerical solutions of periodic orbits attract more and more attention [21–23]. Compared with analytic method, the obtainment of numerical solutions is easier. In this study, the out-plane motion of the spacecraft is ignored, only the periodic motion in the plane of the two primary bodies is considered, employing the numerical method. The planar Lyapunov orbit are obtained by solving the first two equations of Eq. (3) as a Two-Point Boundary Value Problem (TPBVP) [24]. Considering the symmetry of the periodic orbit about axis l1, the problem can be transformed into an Initial Value Problem (IVP) [25]. The initial position is usually chosen on axis l1, thus the initial parameters for Lyapunov orbit are
around the expected orientation can be obtained by adding event constraint conditions. Large amount of efforts can be saved on finding the initial value. Finally, a large number of planning work are carried out to verify the effectiveness of the proposed method. 2. Model The finite-sized rigid body spacecraft is considered in the present paper. In Fig. 1, P1 is the more massive attractor with mass m1, and P2 is the second attractor with mass m2. The mass ratio between P1 and P2 is µ = m2 /(m1 + m2 ) . The distance between P1 and P2 is denoted as d. The dynamics of the spacecraft are modeled under the assumption of planar circular restricted three-body problem (CR3BP), such that: 1) the size and the mass of the spacecraft can be ignored compared with the twoprimary bodies; and 2) P1 and P2 move on a circular orbit at an angular velocity Ω around O, which is the center of mass of P1 and P2. Before specifically developing the work, it is necessary to supply a general framework. Let a primary frame l, rotating with the primaries P1 and P2, be fixed at O. The three unit component vectors are denoted as l1, l2, and l3. l1 is aligned from P1 to P2, l3 is normal to the orbital plane of the primary bodies, and l2 completes a right-handed set of orthogonal axes with l1 and l3. Define the local frame b of the spacecraft to be fixed in the rigid body and attached to the center of mass o. The three unit component vectors denoted as b1, b2, and b3 are along the principal axes of inertia of the spacecraft. The two frames are related by yaw angle , pitch angle , and roll angle , with a rotation sequence of 3-2-1 from the primary frame to the local frame. For a finite-sized rigid spacecraft, a complete description of motion includes the orbital and attitude dynamics. Because of the significantly different scales of variables involved in these two kinds of motion, computational issues may arise. To overcome this problem, the units of mass, length and time are chosen as the sum of mass, the distance and the angular velocity of the two-primary system, respectively, such that m1+m2=1, d=1, and Ω =1. Obviously, the gravitational constant then has the value of unity.
[x y x y ]T = [x 0 0 0 y0 ]T
2.2. The attitude motion Based on the assumption that the spacecraft is influenced only by P1
Denote the position of the spacecraft as [x, y, z] in the primary frame, which is a position vector with respect to the barycenter. Then vectors pointing to the spacecraft from P1 and P2 can be expressed as
r1 = (x + µ ) l1 + yl2 + zl3 1) l1 + yl2 + zl3
(1)
Parameter µ in Eq. (1) is the distance from P1 to O, which is actually equal to the mass ratio after non-dimension. The unit coefficients of vectors r1 and r2 can be expressed as
r1 = r2 =
1 l1 2 l1
+ µ1 l2 + + µ 2 l2 +
1 l3 2 l3
(4)
For a given x0, the orbit is obtained by integrating Eq. (3) until the final y and integral time are changeless within the given tolerance. Fig. 2 shows two L1 planar Lyapunov orbits in the Earth-Moon system. The orbits in Fig. 2 span an amplitude Ay (maximum coordinate value in the l2 direction) of 0.2862 and 0.3331 after non-dimension, respectively.
2.1. The orbital motion
r2 = (x + µ
(3)
Fig. 2. Two planar Lyapunov orbits.
(2) 169
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 3. Phase trajectories of L1 Lagrangian point in the Earth-Moon three-body system.
and P2, while other environmental torques are treated as perturbations. For a rigid spacecraft subjected to gravity gradient torques, the equation of attitude motion of in the CR3BP is described by Euler's equation as
I
+
×I
=
single first-order equation by substituting Eq. (1) and Eq. (6) into Eq. (5)
1¨ =
(5)
g
Here, I is the spacecraft's inertia matrix in diagonal form, denoted as I = diag [I1, I2, I3]. denotes the body frame's angular velocity with respect to the inertial space. The gravity gradient torque g due to P1 and P2 is given as g
=
3m1 3m r1 × Ir1 + 5 2 r2 × Ir2 r15 r2
A=
sin cos 0
0 0 1
2 i=1
3mi ( ri3
2 2 i
Here the coefficient
=
I2
µi2 )sin 2
+ i=1
3m i ri3
i µi
cos 2
(10)
represents the inertia ratio
I1 I3
(11)
The subsequent analysis is carried out based on the attitude equation (10) and the orbital equation (3).
(6)
Eq. (3) and Eq. (5) are solved together to obtain the orbit-attitude solutions of the spacecraft. In this study, rotations and displacements outside the plane P1 and P2 are ignored, that means, z = 0 , 1 = 2 = 0 , and the rotation matrix A is a function of the pitch angle
cos sin 0
1 2
3. Methodology In this section, the attitude dynamic characteristics in phase plane are firstly analyzed. Based on it, boundary conditions, convergent conditions and path constraints are deduced. In the final subsection, the Pseudospectral method for path planning is described.
(7)
Additionally,
l1 = [ cos l2 = [ sin
sin cos
l3 = [0 0 1]T
3.1. Dynamic analysis on phase plane
0]T 0]T
To get the phase trajectories, integrate Eq. (10) over time. Then the following equation can be obtained
(8)
The angular velocity in the body frame with respect to the inertial frame is
=
+
2
+
here,
(9)
Consequently, the description of attitude behavior reduces to a 170
sin(2 + ) = E
(12)
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 4. The Phase Angle
=
and d0 in the Earth-Moon Space.
curve in Fig. 3(b) is greater than that in Fig. 3(a), and it is the same for Fig. 3(c) and (d). According to Eq. (12) and Eq. (13), the amplitude of the critical curve can be expressed as
c12 + c22 2
1 2
c1 =
i=1 2
c2 = i=1
sin =
c1
3mi ( ri3
3mi ri3 , cos
2 i
µi2 )
d0 =
i µi
=
c2
(13)
=
E
(
)Stablemotion around the equilibrium
Critical stable(sinecurvespassingthrough , )divergent motion(wavycurves)
4
(17)
2
0. 2192
(18)
These two solutions correspond to the center of the two special regions in Fig. 4, respectively. As the amplitude d0 approaches zero, the stable regions reduce to a narrow ellipse and finally turn out to be an equilibrium point. Then the closed curves vanish, without any lines passing through the axis = 0 . Thus, all the attitude state including the equilibrium points become unstable, and a diverging pitch motion occurs. Actually, substitute = 0.2 and = 0 into Eq. (12), the following figures can be drawn. Fig. 5(a) shows the phase trajectories near the center of the special region. The small amplitude of the critical curve in Fig. 5(a) limits the stable region to a narrow one. Fig. 5(b) shows the central phase trajectories of the special region. It can be noted that the stable region vanished. Actually, the red lines in Fig. 5(b) correspond to infinite discrete points on axis = 0 . The dynamic characteristics of these points are all balanced. Different from the equilibrium points in Fig. 3, the phase trajectories around these points turn into straight lines instead of closed curves, sine curves or wavy curves. That means the
point(closed curve) E0 = E (
or
x2 = 0. 9373, y2 =
Thestableequlibriumpoint ,
2
x1 = 0. 9373, y1 = 0. 2192
To sum up, the value range of E and the corresponding characteristics on the phase plane are as follows.
Emin =
3 4
d0 determines the size of stable region. The bigger d0 is, the larger the stable region is. Fig. 4 shows the distribution of and d0 in the EarthMoon space. Two regions emerge on the upside and downside of the moon in both Fig. 4(a) and (b). The phase angle changes sharply from −90° to 90° and the amplitude of the critical curve d0 reaches zero in these two regions. Actually, when solving the equation = 0 , two solutions about x and y can be obtained for the Earth-Moon system.
(14)
sin(2 + )) =
(16)
2
If the position of the spacecraft is fixed, E0 is simply decided by . Thus, a larger value of produces a larger amplitude of sine curve. Generally speaking, for a fixed position in the three-body space, four equilibrium points of the pitch angle motion can be found. Their characteristics of stability is decided by the sign of the inertia ratio , and the size of the stable regions is determined by the value of , which ranges from −1 to +1. The above research obtains the characteristics of phase trajectories at a fix position. When it comes to the whole three-body space, it is important to have a general understand of the distribution of equilibrium angle and the size of stable region. According to Eq. (12) and Eq. (16), the phase angle determines the position of equilibrium angle. The equilibrium angle can be expressed as
Note that, E is determined by the position and attitude of the spacecraft, ranging from to positive infinity. For a fix position in the three-body plane, the value of is constant. Then families of phase path can be drawn by changing the value of E. Fig. 3 shows the phase trajectories at L1 Lagrangian point for different inertia ratios. Four equilibrium points can be noted on the phase plane in Fig. 3, with a difference of 90° between neighbors, denoted by red points. The shadow regions represent stable attitude motions. When the attitude state of the spacecraft falls into these regions, it will move on along the closed curve. Thus, the attitude angle will vibrate periodically around the center equilibrium point. However, when the attitude state is beyond these regions, it will change along the wavy curves. In this case, the attitude angle will diverge. If the attitude state is above the shadow region, the angle velocity is positive, and the attitude angle will grow to positive infinity and vice versa. Obviously, the two saddle points are unstable while the other two points are stable. The stable points correspond to the minimum value of E, while the saddle points correspond to a critical value E0. With the increase of E, divergence occurs on the pitch angle motion after it reaches E0. The phase trajectories no longer form a closed curve, turning into a wavy line. Actually, when E is greater than E0, cannot reach the zero point, which means that E is sin(2 + ) . Obviously, the crilarger than the maximum value of tical value of E0 is
E0 = max(
2E0 =
axis) (15)
Compared with Fig. 3(a) and (c), Fig. 3(b) and(d), it can be discovered that the stable equilibrium points between them take a difference of π/2, which can be explained by: a) If > 0 , E0 = and = /4 /2 . b) If < 0 , E0 = /2 . According to and = /4 Eq. (13), for a fixed position, is constant. Thus, changing the sign of leads to a movement of π/2 for stable equilibrium points. Then considering Fig. 3(a) and (b), the amplitude of the critical 171
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 5. Phase trajectories in the two special regions.
Fig. 6. Phase trajectory for Ay = 0.2862 ( = 0.6 ).
balance is easily to be disturbed, and becomes divergent quickly. Actually, the balanced points at the center of the special region are much more unstable than the saddle points discussed before. However, almost
all the common periodic orbits, including Lyapunov orbits and Halo orbits around L1 and L2 Lagrangion points, transit through these regions. According to the results of Knutson et al. [15], spacecraft can 172
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 7. Phase trajectory for Ay = 0.3331 ( = 0.6 ).
hardly remain close to its initial orientation on planar Lyapunov orbit except for certain spacecraft topology and orbit amplitude. Simulations were carried out by Knutson et al. [15], and it was concluded that the spacecraft can remain close to its original attitude ( = = 0 ) on orbits Ay = 0.2862. However, when it comes to orbit Ay = 0.3331, a great change takes place in pitch motion, nearly 180°. Figs. 6 and 7 shows the attitude motions on these two orbits over one revolution in the phase plane. Different figures correspond to different moment in one revolution, as shown in the subtitle. The red curves are the phase trajectories of the spacecraft from the beginning to the
current time. The blue points correspond to the current attitude state, that means the horizontal and vertical coordinates of this point represent the current attitude angle and the attitude velocity, respectively. The black curves are the phase trajectories of the current orbital position. Take Fig. 6 (a) as an example, the horizontal and vertical coordinates of the blue point are the attitude state at 0.34 revolution and the phase trajectory corresponds to the orbital position at 0.34 revolution. It is shown in Fig. 6 that the phase trajectory change a lot from 0.4 revolution to 0.6 revolution, which actually corresponds to the two 173
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 8. The planning phase trajectories for Ay = 0.3331 ( = 0.6 ).
through half a revolution, which directly leads to the divergence of the pitch motion. Therefor, when the phase trajectory moves beyond the stable region, divergence is most likely to take place. 3.2. The path planning schedule Equations of Motion. In this study, we are focused mainly on the global orbit-attitude dynamics of the pitch motion in the circular restricted three-body problem. Thus, the attitude equation (10) and the orbital equation (3) derived above are solved together to achieve the orbit-attitude solutions of the spacecraft. Objective Function. It is usually essential for spacecraft to maintain a certain orientation. Thus, an expected pitch angle is selected, denoted as ˜ (t ) . Then the objective is to plan a maneuver path near the expected orientation as close as possible. In other words, the maximum deviation from the expected orientation should be the smallest. Thus, the objective function can be expressed as
Fig. 9. Spacecraft Configuration example.
min J = max(
(t )
~
(19)
(t ) )
Constraints. The constraints contain boundary conditions and path constraints. In some special cases, the event constraints are also employed. In this study, no artificial control is introduced. The attitude motion is decided simply by the initial attitude and position. The path planning work is actually to find the proper initial state and to achieve the greatest satisfaction. So no boundary conditions of the attitude motion are employed here. The boundary conditions of the orbital motion can be expressed as (20)
x (t 0) = x 0 , y (t 0) = 0, x (t 0 ) = 0, y (t0) = y0
To ensure the stability, the path constraints are set to keep the phase trajectory inside the stable region of the local position. Thus, the path constraints are denoted as
Fig. 10. The Pitch Angle Response on L1 Lyapunov family over one revolution. (Earth-Moon system).
E (t ) =
(t ) 2 +
(t )sin(2 (t ) + (t )) < E0 (t ) =
(t )
(21)
Additionally, periodic solutions can be obtained by adding event constraints
special regions in Fig. 4. During this period, the phase angle changes by about 180°and the amplitude of the critical curve varies from smaller than 1 to larger than 10. The same characteristics can be found in Fig. 7. However, the essential difference between Figs. 6 and 7 is the stability of the attitude motion. The motion on orbit Ay = 0.2862 is always within the scope of stable region, as shown by the black curves. On the contrary, the motion in Fig. 7 enters the unstable region after passing
(t 0) =
(t f )
(t 0) =
(t f )
(22)
Here the subscript 0 means the initial time and f represents the final time. 174
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 11. The Pitch Angle Response on L1 Lyapunov family over multi revolutions. (Earth-Moon system).
3.3. Pseudospectral method
the algebraic equations at LGL nodes and the optimal control problem is transformed into a nonlinear programming (NLP). Then, the NLP can be solved by sequential quadratic programming algorithm. In this paper, the GPOPS is used, which is a PSM Matlab tool box design by Rao [26].
The Legendre pseudospectral method is employed here. The states are estimated with Legendre polynomials Pi that is orthogonal on the interval [-1, 1]. The nth degree Legendre polynomial is
1
Pn (t ) =
2nn!
dn 2 (t dt n
1)n
4. Results
(23)
4.1. Method verification
The states are interpolated at pre-selected node points, and the values at these points become the optimal variables. The dynamic equation (9) and (3) and the path constraints are enforced at the Legendre-Gauss-Labatto (LGL) points. The constraints for the entire interval are assumed to be satisfied when they are satisfied as the node points. Define x (t ) = [x (t ) y (t ) x (t ) y (t ) (t ) (t )]T to be the vector of all states for this problem and let x (t ) = f (x (t ), v~ (t )) represent the corresponding right-hand sides of the dynamic Eq. (9) and (3). For any given time interval [t 0, t f ], the transformation t : [ t0, tf ] = [ 1,1] [t 0, t f ] is given by
(t f
t (t ) =
t 0 ) t + (t f + t 0 )
To begin the path planning work in the three-body problem, an example is employed to verify the effectiveness of the proposed method. The divergence motion on orbit Ay = 0.3331 is planned, which is shown in Fig. 7. To make a comparison, the expected pitch motion is chosen as ˜ (t ) = 0 . Fig. 8 shows the planned results. The planned initial attitude angle is (t 0 ) = 7.15 , just with a little deviation from (t 0 ) = 0 which can be ignored. After planning, the phase trajectory locates in the stable regions over the whole revolution. The maximum deviation from the expected pitch angle is only about 13°. Divergence does not occur here, which demonstrates that the proposed planning method is effective for this case. Additionally, the planned phase trajectories seem to be symmetry, which may indicate a better dynamic characteristic.
(24)
2
Eq. (23) can be used to shift the LGL points to the desired interval. Using the LGL points, the state variables are estimated by
4.2. Mapping the pitch angle response
N
x N (t ( t )) =
x (t )
xi
i (t
)
(25)
i=0
where xi are coefficients, satisfies i (t
i (t
) is the Nth Lagrange polynomials that
1, if i = k 0, if i k
)=
This subsection is focused on planning the attitude motion path along different Lyapunov periodic orbits. As the amplitude of periodic orbits increases, the fully nonlinear dynamics over one revolution become more complicated, and the phase angle and d0 also change greatly. As discussed above, Knutson et al. [15] have carried out a large amount of simulations to investigate the natural attitude response along these orbits. The initial condition was set to be (t ) = (t ) = 0 by them. In the present work, the L1 family is selected again, two orbits of which are displayed in Fig. 2. The planned results will be shown in a map to compare with the natural response by Knutson et al. The variables in the map are chosen as the amplitude of each periodic orbit and the configuration of the spacecraft. The configuration of the spacecraft is described by the inertia ratio , as defined by Eq. (11). Take the model in Fig. 9 as an example to illustrate the relation between the configurations and the inertia ratio . The spacecraft is regarded as an approximate cuboid, the three body vectors are along the principal axes of inertia of the spacecraft and the three side lengths are denoted as a, b and c, respectively. Then Eq. (29) can be obtained. It is observed that when ignoring the rotations and displacements outside of the orbital plane of P1 and P2, the spread of the spacecraft in this direction makes no effects on the attitude motion. Instead of setting initial states, the expected attitude motion is
(26)
The derivatives of the state variables at the LGL nodes can be computed by
x N (t ( t )) =
x (t ( t ))
N
2 (t f
t0 )
xi i=0
i (t
)=
N
2 (t f
t0 )
Dki x i i=0
(27)
The (N + 1) × (N + 1) differential matrix D associated with the Legendre polynomials is PN (tk ) PN (tk )(t f ti )
Dki =
N (N + 1) 4 N (N + 1) 4
0
if k
i
if k = i=0 if k = i=N otherwise
(28)
Thus, the differential equation of the control model is estimated by 175
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 12. Multi-Period Simulation results. (Ay = 0.2037, κ = 0.1).
selected to be ˜ (t ) = 0 for the sake of comparison. Maintaining this expected attitude, different simulations are completed by varying the amplitude of the orbit and the inertia ratio . Fig. 10 is a map of the global deviation from the expected attitude over one revolution. Every point in the map corresponds to a planning work over one revolution, then the maximum bias of the pitch angle from the expected motion is recorded. Different colors in Fig. 10 represent the value of the bias. The computation time for one planning task is averagely 20s with one processor (8 GB RAM).
If = 1, a=0 If 1 < < 0, a < b If = 0, a=b a>b If 0 < < 1, b=0 If = 1,
(29)
On the whole, the planned pitch angle over one revolution on the L1 Lyapunov family maintains closer to the expected orientation than the natural motion with the expected attitude being its initial condition. 176
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 13. Multi-Period Simulation results. (Ay = 0.4556, κ = −0.6).
Besides, all the phase trajectories of the planned results fall into the stable regions on the phase plane. Although the pitch angle grows larger with the increase of Ay, bias larger than 70 deg has not been found within the scope of this study. Compared with the natural motion in Ref. [13], the proposed planning method works well for the fully nonlinear attitude-orbit coupled system. Additionally, a dark horizontal band emerges from the map for a narrow set of Lyapunov orbits, which is consistent with the feature discovered by Knutson. Fig. 11 gives the maps for 2 revolutions and 4 revolutions, respectively. As is mentioned above, the computation time for a one-
revolution planning work costs about 20s. When it comes to multiperiod planning work, a tremendous growth in planning time occurs since the LGL points and the Legendre polynomials increase with the periodicity. The overall characteristics of distribution in Fig. 11 are similar to that in Fig. 10. When the spacecraft is axisymmetric about b3, it can always maintain = = 0 . With the growth of amplitude of the orbit, the pitch angle diverges more greatly from the expected orientation, especially for 4-revolution planning results. The largest angle for 4-revolution reaches 90° in Fig. 11(b). It is worth noting that the dark horizontal band turns less obvious as the planning period 177
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 14. The phase trajectories of periodic solutions on orbit Ay = 0.1817. ( ˜ (t ) = 0 ).
Fig. 15. The phase trajectories of periodic solutions on orbit Ay = 0.3044. ( ˜ (t ) = 0 ).
Fig. 16. The phase trajectories of periodic solutions on orbit Ay = 0.1817. ( ˜ (t ) = /2 ).
Fig. 17. The phase trajectories of periodic solutions on orbit Ay = 0.3044. ( ˜ (t ) = /2 ).
Fig. 18. The phase trajectories of periodic solutions on orbit Ay = 0.1817. ( ˜ (t ) =
178
/2 ).
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo
Fig. 19. The phase trajectories of periodic solutions on orbit Ay = 0.3044. ( ˜ (t ) =
increases. It can be predicted that this band will disappear with the increase of periodicity. Since the multi-period planning work normally takes much time, multi-period simulations are carried out with initial conditions determined by one-revolution planned results, aiming to see whether the one-period planned results work well for multi periods. Two couples of orbits and spacecraft configurations are selected, one of which falls into the regions less than 10° in Fig. 10, and the other one is in the regions more than 50°. The simulation results are shown in Figs. 12 and 13. The multi-period simulations of the selected orbits show good characteristics on maintaining the attitude motion around the expected orientation. In Fig. 12, the phase trajectories keep a circumferential motion in an adjacent area. The pitch angle shows a waving curve less than 10°. Actually, this stable motion can keep going for a much longer simulation time. In Fig. 13, the phase trajectories show a similar-saddle shape. Within the time frame of the study (4 revolution), no divergence occurs. Although the selected orbits cannot represent all the situations, the superiority of the proposed planning method is outstanding when compared with the none-planning simulations.
/2 ).
feature in these phase trajectories is that they are all closed curves, which guarantees the periodicity. Compared with other methods that aim to find periodic solutions, a large amount of efforts on determining the initial guess can be saved by bringing in the expected orientation. 5. Conclusion The path of attitude motion for a rigid body in fully nonlinear Lyapunov periodic orbits is planned in the background of restricted three-body problem. The attitude dynamic characteristics in phase plane are deeply analyzed, and the path constraint conditions are proposed to avoid the divergence of phase trajectories. Then, a path planning method in three-body problem is proposed by employing the Pseudospectral method. Based on this, solving strategy of attitude-orbit coupled periodic solutions is obtained by adding event constraint conditions. A large number of simulations are carried out to explore the attitude behavior across the planar Lyapunov families by changing the configuration of the spacecraft. Results are summarized in maps, which address the question that how the translational motions and the configuration of the satellite affect the planned results. Finally, multiperiod results and attitude-orbit periodic solutions are conducted to make comparison with previous researches. It is shown that the proposed method effectively avoids the divergence of attitude motion, and the bias of pitch angle can be significantly reduced. Besides, a large amount of efforts on determining the initial guess can be saved.
4.3. Periodic solutions As is mentioned before, periodic attitude-orbit solutions can be obtained by adding event constraints described in Eq. (22). It is more complicated that for a selected reference orbit, there may be several periodic solutions, corresponding to different initial pitch angles, which have been discovered by Dayung and Rodney [18]. In their study, not only one-period attitude-orbit solutions but also several multi-period solutions come out for one orbit. In the present study, the expected orientation is also brought in, aiming to select the most satisfied periodic solutions. Based on the existing planning method, event constraints are introduced to guarantee that the pitch angle and angular velocity at the initial time are equal to that at the final time in a revolution. Thus, periodic solutions can be obtained by the planning work. It is worth noting that the attitude-orbit periodic solutions raise a high requirement on the accuracy of periodic orbit. The pre-set tolerance of the iteration to get Lyapunov periodic orbits is 10−5. Two orbits are chosen with an amplitude Ay of 0.1817 and 0.3044 to show the planning periodic solutions in this paper. Three cases are planned with different expected orientations. Case 1: ˜ (t ) = 0 Case 2: ˜ (t ) = /2 /2 Case 3: ˜ (t ) = Figs. 14–19 show the phase trajectories of the planned periodic solutions. Three expected orientations are chosen as ˜ (t ) = 0 , ˜ (t ) = /2 and ˜ (t ) = /2 , corresponding to Case 1, Case 2 and Case 3, respectively. It can be found that the phase trajectories show great difference in shape and size for different orbits, inertia ratios and expected orientations. Generally speaking, the planned trajectories of pitch angle are all moving around the given expected orientation. Certainly, the periodic solutions on one orbit is limited. Not every given expected orientation has a unique corresponding curve. The common
6. Further investigation The planned method in this paper aims to explore the characteristic of attitude motion under the background of CR3BP. But when a certain specific mission is considered, the solutions may be unaccepted. For example, the 70 deg bias over one revolution is not allowed for an observation mission. Thus, the planned method should be improved to consider control [28]. Besides, applying the planning method in the Halo orbit with three dimensions [29] is worthy of further exploration. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 11572345). References [1] R.W. Farquhar, The flight of ISEE-3/ICE: origins, mission history, and a legacy, J. Astronaut. Sci. 49 (1) (2001) 23–74. [2] M.W. Lo, B.G. Williams, W.E. Bollman, D. Han, Y. Hahn, J.L. Bell, B. Barden, Genesis mission design, J. Astronaut. Sci. 49 (1) (2001) 169–184. [3] D.W. Dunham, R.W. Farquhar, G. Gómez, M.W. Lo, J.J. Masdemont (Eds.), “Libration Point Mission, 1978–2002,” Libration Point Orbits and Applications: Proceedings of the Conference, World Scientific Publ. Company, Aiguablava, June 2003. [4] L. Bucci, M. Lavagna, D. Guzzetti, Periodic orbit-attitude solutions along planar orbits in a perturbed circular restricted three-body problem for the Earth-Moon System, Acta Astronaut. 99 (2018) 152–162. [5] D. Guzzetti, K.C. Howell, Attitude dynamics in the circular restricted three-body problem, Astrodynamics 2 (No. 2) (2018) 87–119.
179
Acta Astronautica 163 (2019) 168–180
M.-y. Yu and Y.-z. Luo [6] T.R. Kane, E.L. Marsh, Attitude stability of a symmetric satellite at the equilibrium points in the restricted three-body problem, Celest. Mech. Dyn. Astron. 4 (No. 1) (1971) 78–90. [7] T.R. Kane, P.M. Barba, Attitude stability of a spinning satellite in an elliptic orbit, J. Appl. Mech. 33 (No. 2) (1966) 402–405. [8] W.J. Robinson, Attitude stability of a rigid body placed at an equilibrium point in the restricted problem of three bodies, Celest. Mech. Dyn. Astron. 10 (No. 1) (1974) 17–33. [9] W.J. Robinson, The Restricted problem of three bodies with rigid dumb-bell satellite, Celest. Mech. Dyn. Astron. 8 (2) (1973) 323–330. [10] A. Abad, M. Arribas, A. Elipe, On the attitude of a spacecraft near a Lagrangian point, Astron. Inst. Czechoslovakia Bull. 40 (5) (Sept. 1989) 302–307. [11] E. Brucker, P. Gurfil, Analysis of gravity-gradient perturbed rotational dynamics at the collinear Lagrange points, J. Astronaut. Sci. 55 (3) (2007) 271–291. [12] B. Wong, R. Patil, A. Misra, Attitude dynamics of rigid bodies in the vicinity of Lagrangian points, J. Guid. Control Dyn. 31 (1) (2008) 252–256. [13] M. Lara, J. Peláez, C. Bombardelli, F.R. Lucas, M. Sanjurjo-Rivo, D. Curreli, E.C. Lorenzini, D.J. Scheeres, Dynamic stabilization of L2 periodic orbits using attitude-orbit coupling effects, J. Aerosp. Eng. 4 (1) (2012) 73–82. [14] Y.H. Meng, R. Hao, Q.F. Chen, Attitude stability analysis of a dual-spin spacecraft in halo orbits, Acta Astronaut. 99 (2014) 318–329. [15] A.J. Knutson, D. Guzzetti, K.C. Howell, Attitude responses in coupled orbit-attitude dynamical model in Earth–moon Lyapunov orbits, J. Guid. Control Dyn. 38 (7) (2015) 1264–1273. [16] D. Guzzetti, K.C. Howell, “Coupled Orbit-Attitude Dynamics in the Three-Body Problem: A Family of Orbit-Attitude Periodic Solutions,” AIAA/AAS Astrodynamics Specialist Conference, Aug. 2014 AIAA Paper 2014-4100.
[17] D. Guzzetti, K.C. Howell, Natural periodic orbit-attitude behaviors for rigid bodies in three-body periodic orbits, Acta Astronaut. 130 (2017) 97–113. [18] D. Koh, R.L. Anderson, Periodic orbit-attitude solutions in the planar elliptic restricted three-body problem, J. Guid. Control Dyn. 41 (3) (2018) 644–657. [19] H. Poincaré, New Methods of Celestial Mechanics. History of Modern Physics and Astronomy, American Institute of Physics, 1892. [20] D.L. Richardson, Analytical construction of periodic orbits about the collinear points, Celestial Mech. 22 (3) (1980) 241–253. [21] J. Breakwell, J. Brown, The “Halo” family of 3-dimensional periodic orbits in the Earth-Moon restricted 3-body problem, Celestial Mech. 20 (1979) 389–404. [22] M. Xu, T. Tan, S. Xu, Research on the transfers to Halo orbits from the view of invariant manifolds, Sci. China Phys. Mech. Astron. 55 (4) (2012) 671–683. [23] T. Luo, M. Xu, Y.F. Dong, Natural formation flying on quasi-halo orbits in the photogravitational circular restricted three-body problem, Acta Astronaut. 149 (2018) 35–46. [24] P.B. Bailey, L.F. Shampine, P.E. Waltman, Nonlinear Two Point Boundary Value Problems, Academic press, 1968. [25] J.D. Lambert, Numerical Methods for Ordinary Differential Systems: the Initial Value Problem, John Wiley & Sons, Inc., 1991. [26] A.V. Rao, D.A. Benson, Algortihm 902: GPOPS, a Matlab software for solving multiple-phase optiaml control using the gauss pseudospectral method, ACM Trans. Math Software 37 (2) (2010) 1–39. [28] M. Xu, Y.H. Jia, S.J. Xu, Attitude description and dynamics modeling for Halo-orbit spacecraft, J. Beijing Univ. Aeronaut. Astronaut. 33 (10) (2007) 1166–1169. [29] M. Xu, Y.Y. Liang, T. Luo, X. Sun, Attitude pointing schemes and spacecraft configurations for libration-point-orbit spacecraft, Aero. Sci. Technol. 65 (2017) 1–8.
180