Next Article in Journal
Architecture of Distributed Control System for Gearbox-Free More Electric Turbofan Engine
Next Article in Special Issue
Free-Vertex Tetrahedral Finite-Element Representation and Its Use for Estimating Density Distribution of Irregularly-Shaped Asteroids
Previous Article in Journal
Numerical Investigation of a Rectangular Jet Exhausting over a Flat Plate with Periodic Surface Deformations at the Trailing Edge
Previous Article in Special Issue
Two-Stage Pursuit Strategy for Incomplete-Information Impulsive Space Pursuit-Evasion Mission Using Reinforcement Learning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A New 3D Shaping Method for Low-Thrust Trajectories between Non-Intersect Orbits

1
School of Aerospace Engineering, Tsinghua University, Beijing 100084, China
2
Beijing Institute of Aerospace Systems Engineering, Beijing 100076, China
*
Author to whom correspondence should be addressed.
Submission received: 1 September 2021 / Revised: 8 October 2021 / Accepted: 11 October 2021 / Published: 24 October 2021
(This article belongs to the Special Issue Spacecraft Dynamics and Control)

Abstract

:
This paper proposes a new shape-based method in spherical coordinates to solve three-dimensional rendezvous problems. Compared with the existing shape-based methods, the proposed method does not need parameter optimization. Moreover, it improves the flexibility of orbit fitting, greatly reduces the velocity increment and maximum thrust acceleration, and ensures the orbit safety to a certain extent. The shaping function can provide the initial estimate for numerical trajectory optimization and improve the convergence rate in a certain range when combined with the normalization method. The superiority of the proposed method over the existing methods is demonstrated by two numerical examples. Its effectiveness at initial estimation generation in the indirect optimization of a low-thrust trajectory is demonstrated by the third example.

1. Introduction

The continuous low-thrust propulsion system has attracted great attention because of its higher specific impulse and consequently less fuel consumption compared with impulsive propulsion systems. There are three main categories of optimization methods to optimize the low-thrust trajectory: the direct methods [1,2,3], indirect methods [3,4,5,6], and hybrid methods [7,8,9]. For the gradient-based versions of these methods, a good initial value estimate is key to improving the convergence of optimization. To provide the initial guesses, especially for a general three-dimensional (3D) rendezvous problem, the shape-based method can be used to approximate the near-optimal trajectory through shaping function design. In particular, when the shape is used as a nominal trajectory and then combined with the indirect method, it can greatly improve the computational efficiency [10].
At present, the shape-based methods can be divided into two categories according to whether their shape functions consist of free optimization parameters: one has no parameter optimization, and the other has to resort to the help from parameter optimization. The former has obvious advantages in terms of computational speed and efficiency. Petropoulos and Longuski [11] proposed the first shape-based method without optimization, which uses an exponential sinusoid to describe a coplanar trajectory. It is simple in form but cannot meet the boundary conditions of position and velocity at the same time, so it cannot solve rendezvous problems. Then, the inverse polynomial method was introduced by Wall and Conway [12]. This has enough coefficients to satisfy the in-plane boundary conditions. Because of its limited number of extrema, it cannot solve multirevolution problems. Accordingly, Wall and Pols [13] used the cosine inverse polynomial (cosine IP) to design the coplanar trajectory. In order to solve the problem of 3D transfer and rendezvous problems, Wall [14] and Novak [15] both designed the out-of-plane component; however, this only performs well when the elevation angle is small. Later, Xie [16,17] and Zeng [18] proposed methods that are suitable for transfer and rendezvous problems with any elevation angle, respectively. The difference between these two methods is that Xie’s method [17] designs the in-plane and out-of-plane shapes separately and selects the polynomial as the shaping function, while Zeng’s method [18] does not design the shape separately in three dimensions and uses a Fourier series as the shaping function. Both of these approaches can provide solutions for transfer and rendezvous problems with large elevation angles. However, the velocity increment Δ V and maximum thrust acceleration T a , m a x of the two methods are much larger than those of the 3D inverse polynomial method, which makes the trajectory designed by these two methods far from the optimal trajectory.
There are two main contributions of this paper. Firstly, a new shape-based 3D method is proposed for rendezvous problems, which considers the advantages of Novak’s method [15] and Zeng’s method [18] in the design and improves some deficiencies existing in these two methods. It widens the application scope of the shape-based method and reduces the fuel consumption. The proposed method satisfies the boundary conditions and the equation of motion constraints and does not require parameter optimization. It is able to effectively reduce the velocity increment Δ V and maximum thrust acceleration T a , m a x and provide good initial guesses for the local optimization of low-thrust trajectories. Secondly, adjoint normalization [19] is added into the process of combining the shape-based method with the indirect method, which greatly improves the convergence rate.
This paper is organized as follows. In Section 2, a new shape-based method is proposed. Firstly, the dynamic model of a spacecraft and the 3D trajectory models in spherical coordinates are given, and then the shaping functions are determined. Section 3 introduces the method of applying adjoint normalization into the process of combining the shape-based method with the indirect method in cylindrical coordinates. The specific examples of rendezvous problems and the shape-based method in combination with optimization problems are described in Section 4. Finally, the conclusion is deduced in Section 5.

2. Shape-Based Method for 3D Orbital Rendezvous

In this section, a shape-based method that approximates the low-thrust trajectory with the shaping functions parameterized in the spherical coordinate frame is introduced for the three-dimensional (3D) orbital rendezvous. The boundary conditions and dynamical constraints of the rendezvous are automatically satisfied by developing a new shape-based approach based on the cosine IP method [13] and Zeng’s method [18] empirically to reduce fuel consumption.

2.1. Low-Thrust Dynamics and Shape Design in Spherical Coordinates

The spacecraft in an interplanetary rendezvous is subject to the central gravitational force of the Sun and propelled by a low-thrust propulsion system. For the 3D trajectory design, the corresponding dynamics are given [17] as follows:
r ¨ = μ r r 3 + u
where r = r cos φ cos θ , r cos φ sin θ , r sin φ T is the position vector of a spacecraft with r = r in the spherical inertial frame, u is the low-thrust acceleration, and μ is the gravitational parameter of the Sun. The spherical coordinate frame T O ; r , φ , θ is defined as Figure 1, where r is the distance from the Sun to the spacecraft, θ is the azimuthal angle, and φ is the elevation angle. Equation (1) can also be written as
θ ˙ 2 r + θ ¨ r = μ r r 3 + u
where * and * denote the derivative and double derivative with respect to θ , respectively. The trajectory parameterized with r θ , φ θ , t θ will be given by the shaping functions, and the low-thrust acceleration can then be obtained by Equation (2) to approximate the rendezvous.
The shape designed by Novak [15] empirically has a smaller Δ V and T a , m a x in most coplanar transfer and rendezvous cases compared with other methods. However, for the cases with a large elevation angle φ , it may lead to t < 0 or φ > π / 2 and thus cannot obtain a reasonable trajectory. Inspired by the shaping functions proposed in [18], an improved interpolation between the initial and target orbits is proposed here to design the elevation angle φ θ . Combined with the r θ designed in [15] by introducing a mid-plane, the proposed shape can broaden the application range to the 3D rendezvous trajectory with a large elevation angle and reduce the velocity increment and maximum low-thrust acceleration to a certain extent empirically.

2.1.1. Shaping the Elevation Angle

The function φ θ is defined according to the direction of spacecraft position e r = e r , x , e r , y , e r , z T , viz.,
φ θ = arctan e r , z e r , x 2 + e r , y 2 π 2 , π 2
where the direction vector e r (for the convenience of derivative calculation, e r is not a unit vector) is obtained by an interpolation between the unit vectors in the initial and target orbital planes (i.e., e r 1 and e r 2 , respectively) in combination with a shaping function ϕ θ :
e r = ϕ e r 1 + ( 1 ϕ ) e r 2
This shaping function determines the elevation angle and is given empirically below. Note that the unit vectors e r j , j = 1 , 2 represent the directions at θ constrained in the initial and target orbital planes, respectively, which can be expressed as
e r j = cos f j p j + sin f j q j
where f j is the true anomaly at θ , p j is directed along the Laplace–Runge–Lenz vector, and q j is along the cross product of vector of angular momentum and p j . The classic orbital elements of the initial and final orbits are known in advance, and the direction vectors p j and q j are given by
p j = p j x p j y p j z = cos ω j cos Ω j sin ω j sin Ω j cos i j cos ω j sin Ω j + sin ω j cos Ω j cos i j sin ω j sin i j q j = q j x q j y q j z = sin ω j cos Ω j cos ω j sin Ω j cos i j sin ω j sin Ω j + cos ω j cos Ω j cos i j cos ω j sin i j
where i j is the inclination, Ω j is the right ascension of the ascending node, and ω j is the argument of the perigee of tje corresponding orbital plane. There is a projection relationship between f j and the variable θ in the spherical coordinate system; i.e., tan θ = e r , y / e r , x . Thus, f j can be expressed as a function of θ :
f j = atan 2 p j , x sin θ p j , y cos θ , q j , y cos θ q j , x sin θ π , π

2.1.2. Shaping the Radius

The cosine IP method [13] is used to design the two-dimensional (2D) trajectory due to its advantages of reduced fuel consumption and a simple form in practical experience. Therefore, the radius of the spacecraft in the spherical coordinates is shaped by the following function:
r = 1 k 0 + k 1 θ + k 2 θ 2 + ( k 3 + k 4 θ ) cos θ + ( k 5 + k 6 θ ) sin θ
where the coefficients k j and j = 0 , 1 , , 6 are unknown free parameters used to satisfy the boundary conditions and the constraint on the time of flight.
Equation (8) is firstly proposed for the case of a rev2D trajectory, where θ is the parameter in polar coordinates. For the 3D case with a large elevation angle, the corresponding true anomaly shown in Figure 2 changes sharply with respect to the azimuthal angle θ around the descending/ascending nodes. In this case, the radius obtained by Equation (8), which is intentionally designed for the 2D case, is unreasonable around these nodes. Thus, the trajectory is preferred to be parametrized in terms of the true anomaly in the transfer orbital plane compared with the azimuthal angle in the x o y plane. However, considering the difficulty in determining the relationship between the true anomaly of the spacecraft and the azimuthal angle, a compromise strategy is proposed here. Instead of the time-varying transfer orbital plane, a fixed middle orbital plane between the initial and target orbits is used as a reference plane, which can ensure that the angle between the reference plane and the x o y plane is within 45 and avoid the change of the curvature direction of the shape as much as possible when the elevation angle is large. The corresponding true anomaly is given by
f = atan 2 p x sin θ p y cos θ , q y cos θ q x sin θ
As shown in Figure 2, the unit vectors p and q point to the semi-major and semi-minor axes of the middle orbital plane, respectively, which are obtained by
p = p x p y p z = cos ω 1 + ω 2 2 cos Ω 1 + Ω 2 2 sin ω 1 + ω 2 2 sin Ω 1 + Ω 2 2 cos i 1 + i 2 2 cos ω 1 + ω 2 2 sin Ω 1 + Ω 2 2 + sin ω 1 + ω 2 2 cos Ω 1 + Ω 2 2 cos i 1 + i 2 2 sin ω 1 + ω 2 2 sin i 1 + i 2 2 q = q x q y q z = sin ω 1 + ω 2 2 cos Ω 1 + Ω 2 2 cos ω 1 + ω 2 2 sin Ω 1 + Ω 2 2 cos i 1 + i 2 2 sin ω 1 + ω 2 2 sin Ω 1 + Ω 2 2 + cos ω 1 + ω 2 2 cos Ω 1 + Ω 2 2 cos i 1 + i 2 2 cos ω 1 + ω 2 2 sin i 1 + i 2 2
Accordingly, Equation (8) becomes
r = 1 k 0 + k 1 f + k 2 f 2 + ( k 3 + k 4 f ) cos f + ( k 5 + k 6 f ) sin f
By substituting Equation (9) into Equation (11), the shaping function r θ is finally formulated for the 3D rendezvous trajectory.

2.1.3. Shaping the Time

The unit vectors e r , e o , e h define the radial–orthoradial–out-of-plane coordinate system [15]. We project the dynamic equations to the tangential–normal–out-of-plane reference frame, which is defined by the unit vectors e t , e n , e h , and then Equation (2) becomes
u = | u t u n u h = μ r 2 e r · e t + θ ¨ v ˜ · e t + θ ˙ 2 a ˜ · e t μ r 2 e r · e n + θ ˙ 2 a ˜ · e n θ ˙ 2 a ˜ · e h
where v ˜ and a ˜ are the velocity and acceleration vectors in the Cartesian coordinates, respectively, and can be written as
v ˜ = d r d θ = P S C v ˜ S , a ˜ = d 2 r d θ 2 = P S C a ˜ S
v ˜ S = r r cos φ r φ , a ˜ S = r r φ 2 + cos 2 φ 2 r cos φ 2 r φ sin φ 2 r ϕ + r φ + sin φ cos φ
P S C = cos θ cos φ sin θ cos θ sin φ sin θ cos φ cos θ sin θ sin φ sin φ 0 cos φ
where P S C is the coordinate transformation matrix from the spherical coordinates to the Cartesian coordinates, and v ˜ S and a ˜ S denote the velocity and acceleration vectors in spherical coordinates, respectively. The expressions of r , r , φ , and φ are obtained by the chain rule:
r = d r d f f
r = d 2 r d f 2 f 2 + d r d f f
φ = φ ϕ ϕ + φ f 1 f 1 + φ f 2 f 2
φ = 2 φ ϕ 2 ϕ 2 + 2 φ f 1 2 f 1 2 + 2 φ f 2 2 f 2 2 + 2 2 φ f 1 f 2 f 1 f 2 + 2 φ f 2 ϕ f 2 ϕ + 2 φ f 1 ϕ f 1 ϕ + φ ϕ ϕ + φ f 1 f 1 + φ f 2 f 2
The detailed formulation of this is given in the Appendix.
Continuous tangent thrust has the advantages of reducing energy consumption in orbit transfer and avoiding the need to change the thrust direction in the orbit plane. It is also convenient for time shaping.Therefore, assuming that there is no normal component of u , which is equivalent to u n = 0 , from Equation (12), one can obtain
t = 1 θ ˙ = D r 2 μ
the expression of D is
D = r + 2 r 2 r + r φ φ sin φ cos φ φ 2 + cos 2 φ + r φ 2 + cos 2 φ
The sign of D can represent the curvature direction of the trajectory. Equation. (20) requires that D must be positive, which means that the trajectory must bend towards the central celestial body.
Combining Equation (12) and Equation (20), the velocity increment during the orbit transfer can be obtained as
Δ V = θ 1 θ 2 u ( θ ) θ ˙ d θ
The angular acceleration θ ¨ is contained in u ( θ ) , which can be calculated as
θ ¨ = θ ˙ 3 D r 2 + 2 D r r 2 μ D r 2
where
D = r + 4 r r r 2 r 3 r 2 + r φ + r φ φ sin φ cos φ φ 2 + cos 2 φ + r φ φ φ cos 2 φ + φ sin 2 φ φ 2 + cos 2 φ + 2 φ φ 2 φ cos φ sin φ r r φ φ sin φ cos φ φ 2 + cos 2 φ 2 + r φ 2 + cos 2 φ
r and φ are also obtained by the chain rule, the details of which are shown in Appendix A.

2.2. Determination of Shaping Functions

To determine the interpolation function ϕ and the seven unknown coefficients for the radius shape, the trajectory boundary constraints and non-intersect constraint are analyzed in this subsection. The boundary constraints on the rendezvous trajectory consist of 11 shape boundary constraints (including the initial states r 1 , φ 1 , r 1 , φ 1 , t 1 T , the final states r 2 , φ 2 , r 2 , φ 2 , t 2 T , and the transfer time Δ t = ( t 2 t 1 ) and the fixed initial and final independent variables θ 1 , θ 2 . Meanwhile, the unknown coefficients for the radius shape are designed to satisfy the constraints on r 1 , r 1 , t 1 , r 2 , r 2 , t 2 , Δ t , whereas the constraints on φ 1 , φ 1 , φ 2 , φ 2 will be satisfied by the interpolation function in combination with another four free parameters.

2.2.1. Determination of Interpolation Function

The constraints on the elevation angle are satisfied by
ϕ θ 1 = 1 , ϕ θ 2 = 0 , ϕ θ θ = θ 1 = 0 , ϕ θ θ = θ 2 = 0
The interpolation function takes the value of 1 when θ = θ 1 and 0 when θ = θ 2 . Since the initial and final velocities of the spacecraft out of orbital planes are zeros, the derivatives ϕ are zeros at θ 1 and θ 2 .
Suppose that the initial and target orbits do not intersect for any θ θ 1 , θ 2 ; the ϕ angle of the transfer orbits is designed between that of the initial and target orbits (viz. trajectory safety in [18]), which can avoid collision and save energy to a certain extent. Therefore, a monotonic decreasing function ( ϕ < 0 when θ θ 1 , θ 2 ) is designed, and its evolution direction is from the initial plane to the target plane.
There are many forms of ϕ that can satisfy the above constraints. Considering a reduction of the total velocity increment as much as possible, the transfer orbits have to raise the elevation angle as much as possible at the place with the largest radius and not change the elevation angle as far as possible in the place with a small radius. Therefore, the interpolation functions are given empirically for the cases r 1 < r 2 and r 1 > r 2 , respectively:
If r 1 < r 2 ,
ϕ = a + b β + c β n 1 + d β n 2 n 1 n 2 n 1 , n 2 > 1 ;
else if r 1 > r 2 ,
ϕ = a + b β + c ( β + 1 ) n 3 + d ( β + 1 ) n 4 n 3 n 4 n 3 , n 4 < 1
where β = θ θ 0 θ f θ 0 [ 0 , 1 ] is the normalized parameter, and r 1 and r 2 are the radial distances of the departure and arrival positions, respectively.
Since the coefficient matrix of a , b , c , and d is linear, they can be obtained analytically by satisfying the boundary conditions:
if r 1 < r 2
d = n 1 n 2 n 1 c = n 2 n 2 n 1 b = 0 a = 1
Similarly, if r 1 > r 2 , substituting Equation (25) into Equation (27) gives
d = 1 1 2 n 4 + n 4 + 2 + 2 n 4 1 + 2 n 3 n 3 n 4 n 3 2 + 2 n 3 c = 2 n 4 1 1 n 4 2 n 3 1 1 n 3 d b = n 3 c n 4 d a = 1 c d
where n 1 , n 2 , n 3 , and n 4 are the shaping parameters that shape the trajectory. If the semimajor axes, eccentricities, and inclinations of the initial and final orbits are similar, then the absolute values of n 1 , n 2 , n 3 , and n 4 can take small values; for example, from 1 to 10. If the difference is large, then the absolute values of the shaping parameters can take larger values, such as from 10 to 30. Figure 3 shows the evolution of ϕ for these different shaping parameters. In practical problems, the shaping parameters are selected by experience, and usually the transfer orbit is designed to raise the elevation angle as much as possible at the place with the largest radius (apogee) in order to reduce fuel consumption. The values of the coefficients selected by experience already have a good effect in the examples, and the results are not sensitive to these coefficients, so there is no further optimization in this paper.

2.2.2. Determination of Coefficients in Radius

It is mentioned in Section 2.1 that the seven coefficients of Equation (11) are used to satisfy seven constraints.
At the initial point, from Equation (11), it follows that
r θ 1 = r 1 θ 1
r θ 1 = 1 r ˙ 1 t 1
t θ 1 = 1 θ ˙ 1 t 1
Similarly, at the departure point, it yields
r θ 2 = r 2 θ 2
r θ 2 = 1 r ˙ 2 t 2
t θ 2 = 1 θ ˙ 2 t 2
The last constraint is on the time of flight:
θ 1 θ 2 D r 2 μ = t 2 t 1
Using Equations (30)–(36), the coefficients k j , j = 0 , 1 , , 6 can be obtained.
All these seven coefficients are required to be solved numerically. The direct estimation of the initial value of the coefficients will lead to a large amount of calculation and difficult convergence. However, if one adjusts the range of the independent variables from θ 1 , θ 2 to 0 , θ 2 θ 1 , the analytic solutions of some coefficients can be given according to some boundary conditions. Solving Equations (30) and (31) yields
k 0 = 1 r 1 θ 0 k 3 k 1 = r ˙ 1 θ 0 r 1 θ 0 2 k 4 k 5
Then, the number of unknown coefficients will be reduced to 5, which can greatly reduce the computational burden. However, five unknowns and five nonlinear equations correspond to multiple sets of solutions. In this paper, one set of reasonable solutions is sufficient. The fsolve function in Matlab is used to solve k j , j = 2 , 3 , , 6 numerically. We set the initial values of all coefficients to 0, then solve the nonlinear equations, and finally obtain a set of reasonable solutions.

3. Application in Indirect Trajectory Optimization

The combination of the indirect method and parameter continuation method [20] can avoid the estimation of the initial value and ensure the convergence of the optimal control problems, but it is necessary to sacrifice the computational efficiency. Directly estimating the initial value of the adjoint variables can improve the computational efficiency, but it is difficult to realize because the indirect method for the trajectory optimization is sensitive to the initial values of the adjoint variables. An adjoint estimation method is proposed in [21] by linearization around a shape-based path. The estimation method improves the convergence rate of the indirect method for coplanar rendezvous trajectory optimization, where the maximum inclination is only 13 . 54 . For a general 3D rendezvous with large out-of-plane motion, the estimation performance needs further improvement. In this section, the proposed shape-based method is employed to estimate the initial value of the adjoint variables in combination with the estimation method, of which the convergence rate greatly depends on the shaping approximation. The adjoint normalization is also used to improve the convergence rate. Note that there are two essential evaluation indexes for the shape designed: T a , m a x and Δ V . Generally, the smaller the values of T a , m a x and Δ V , the closer the trajectory designed by the shape-based method is to that optimized by the indirect method. The low-thrust trajectory optimization model is first introduced, and then the estimation using the proposed shape in the spherical coordinate system is briefly described.

3.1. Optimization Model

The nonlinear optimal control problem with a fixed flight time is given as follows:
Minimize
J = T m a x I s p g 0 t 0 t f u d t
Subject to
x ˙ ( t ) = A ( x , t ) + T m a x u m B α m ˙ ( t ) = T m a x u I s p g 0 x ( t 0 ) = x 0 , m ( t 0 ) = m 0
where x R 6 is the state vector of the spacecraft, m is the mass of the spacecraft, A ( x , t ) R 6 and B ( x , t ) R 6 × 3 depend on the coordinate system used, u = T / T m a x [ 0 , 1 ] is the engine thrust ratio, α is the unit vector of thrust direction, g 0 is the standard acceleration of gravity at sea level, and I sp is the thruster-specific impulse. The physical meaning of J in Equation (38) is the total fuel consumption. It is difficult to solve the fuel-optimal problem directly because of the bang-bang structure. The homotopy method provides an idea by adding a homotopic parameter ε to the index [19,22]. The new performance index is
J = T m a x I s p g 0 t 0 t f [ u ε u ( 1 u ) ] d t
By reducing ε from 1 to 0, the energy-optimal problem can be connected with the fuel-optimal problem, which is called the continuation technique for homotopic problems. In the process of reducing ε , the solution of the previous step is always taken as the initial estimated value of the current step. Only when ε = 0 will the difficult bang-bang control appear. In other cases, the optimal control law is continuous.
In the homotopic process, when the parameter ε is close to 0, the right hand sides of the ODEs vary rapidly around the switching points. The numerical integrator should be carefully performed when crossing these points. Here, for simplicity of operation, a fixed step integrator with switching function detection rather than an adaptive step size integrator is used [19].
Obviously, when ε = 1, Equation (40) becomes the energy-optimal index:
J = T m a x I s p g 0 t 0 t f u 2 d t
In the next subsection, based on the problem of energy-optimal problem, the combination of the indirect method and shape-based method is discussed.

3.2. Adjoint Estimation Based on Shape-Based Method

This subsection is similar to the problem formulation in [21], which is established accordingly in the cylindrical coordinates [23] and then combined with the indirect method of the optimization problem to obtain the corresponding solution.
The cylindrical coordinates are selected because experience shows that solving the optimization problem in the cylindrical coordinates results in a higher convergence rate than that in the spherical coordinates. Therefore, in practice, it is necessary to convert the shape (nominal trajectory) from the spherical coordinates to the cylindrical coordinates through the coordinate conversion method.
We denote the combination coordinates of position and velocity expressed in the cylindrical coordinates r , θ , z , v r , v θ , v z T by x , and their corresponding adjoint variables λ r , λ θ , λ z , λ v r , λ v θ , λ v z T by λ x . After linearizing the motion equation near the nominal trajectory which is designed in Section 2 and taking the energy-optimal problem as the one that needs to be solved, the motion equations and Euler–Lagrange equations with θ as the independent variable are derived:
d δ x d θ d λ x d θ = F ( θ ) δ x λ x H ( θ )
where
F ( θ ) = M ( x nom ) 0 3 × 3 0 3 × 3 0 3 × 3 I s p g 0 T m a x 2 m 2 I 3 × 3 0 6 × 6 M T ( x nom ) 1 θ ˙ , H ( θ ) = 0 3 × 1 a nom 0 6 × 1 1 θ ˙
where x nom and a nom represent the state and acceleration of the nominal trajectory, respectively. M ( x n o m ) is the Jacobian matrix of the central gravitational acceleration:
M ( x nom ) = 0 0 0 1 0 0 v θ / r 2 0 0 0 1 / r 0 0 0 0 0 0 1 v θ / r 2 + μ 3 r 2 R 2 / R 5 0 3 μ r z / R 5 0 2 v θ / r 0 v r v θ / r 2 0 0 v θ / r v r / r 0 3 μ r z / R 5 0 μ 3 z 2 R 2 / R 5 0 0 0
where R = r 2 + z 2 represents the distance between the central celestial body and the spacecraft.
The solution of Equation (42) is
δ x ( θ f ) λ x ( θ f ) = Φ ( θ f , θ 0 ) δ x ( θ 0 ) λ x ( θ 0 ) θ 0 θ f Φ ( θ f , τ ) H ( τ ) d τ
where Φ ( θ f , θ 0 ) can be derived by integrating the differential equation:
d Φ ( θ , θ 0 ) d θ = F ( θ ) Φ ( θ , θ 0 ) , Φ ( θ 0 , θ 0 ) = I 12 × 12
We denote Φ ( θ f , θ 0 ) as
Φ ( θ f , θ 0 ) = Φ 1 , 1 Φ 1 , 2 Φ 2 , 1 Φ 2 , 2
Let z 1 , z 2 T = θ 0 θ f Φ ( θ f , τ ) H ( τ ) d τ . Then, the initial adjoint variables can be expressed as
λ x ( θ 0 ) = Φ 1 , 2 1 ( δ x ( θ f ) Φ 1 , 1 δ x ( θ 0 ) + z 1 )
The nominal trajectory designed in Section 2 has already met the boundary conditions, which means δ x ( θ 0 ) = 0 and δ x ( θ f ) = 0, so Equation (48) can be simplified as
λ x ( θ 0 ) = Φ 1 , 2 1 z 1
The mass adjoint variable λ m t 0 is given by integrating the following equation:
λ ˙ m = I s p g 0 T m a x 2 m 3 λ v 2

3.3. Adjoint Normalization

In practice, it can be found that the method in Section 3.2 cannot significantly improve the convergence rate when the shapes of the initial and target orbits differ greatly. This subsection introduces the normalization method [19] based on Section 3.2 to improve the convergence rate.
Multiplying the index Equation (41) by a positive multiplier λ 0 does not change the essence of the problem. Consequently, the updated Hamiltonian function is expressed as
H = λ 0 T m a x I s p g 0 u 2 + λ x _ n o r m T M ( x n o m ) δ x + λ v _ norm T T m a x u m α a n o m λ m _ norm T m a x u I s p g 0
where λ v _ n o r m is the group of the last three components of λ x _ n o r m .
Therefore, the normalized adjoint variables λ x and λ m can be expressed as
λ x _ n o r m = λ 0 λ x
λ m _ n o r m = λ 0 λ m
where the positive multiplier λ 0 is obtained according to the normalization condition:
λ 0 = 1 λ x t 0 2 + λ m 2 t 0 + 1
In practice, the initial mass adjoint is randomly selected between 0 and the estimated value λ m _ norm t 0 [21].

4. Numerical Examples

In this section, three orbit rendezvous examples are given to verify the feasibility and advantages of the shape-based method proposed. The Novak method has the limitation of the elevation angle, which means that it fails to obtain a feasible solution when the elevation angle is large. The proposed method is not affected by the inclination of the initial and target orbits and can generate an arbitrary 3D rendezvous trajectory. Compared with the Zeng method, the proposed method can greatly reduce the velocity increment and the T a , m a x , which indicates that the rendezvous trajectory generated by the proposed method is closer to the optimal trajectory. The first two examples illustrate the advantages mentioned above. As for the third example, the trajectory generated by the proposed method is combined with the indirect method and the adjoint normalization to prove that it can broaden the convergence range and improve the convergence efficiency of rendezvous problems. All the examples were performed on a personal desktop with an Intel Core i7-7700 CPU of 3.6 GHz and 8.00 GB of RAM and with MATLAB R2018a.
Nondimensionalization was performed by setting the length unit to the astronomical unit (AU), the time unit, denoted by TU, to 1 / 2 π of the period of a circular orbit with radius 1 AU, and the mass unit, denoted by MU, to the initial mass of the spacecraft. Furthermore, the new velocity unit is denoted by VU.

4.1. Rendezvous Mission A from Inner Orbit to Outer Orbit

The initial and target orbit elements of mission A are listed in Table 1. The elevation angle is set to 60 to verify the feasibility and efficiency of the proposed shape for the 3D rendezvous. The semi-major axis is set from 1 to 4, and the eccentricity is set from 0 to 0.1, which can be used to prove the ability to perform an elliptic rendezvous. Assume that the revolution number N is 3, the fixed flight time t f is 9 years (56.5089 TU), and the shaping parameters are 10 and 20 for n 1 and n 2 , respectively. The shaping parameters are chosen according to experience in order to raise the elevation angle at the place with the largest radius as much as possible.
Then, the trajectories obtained by the three methods are shown by Figure 4 and Figure 5. The velocity increment and the maximal thrust acceleration T a , m a x of each method are listed in Table 2. Compared with Novak’s method and Zeng’s method, the proposed method saves Δ V by 43.22 % and 84.44 % , respectively, which greatly reduces the fuel consumption in rendezvous problems. The maximal thrust acceleration T a , m a x is also significantly reduced from Zeng’s method, at 3.458, and Novak’s method, at 0.2943, to the proposed method’s result of 0.1527.
The longer the rendezvous time, the easier it is to solve the rendezvous trajectory, which can test the flexibility of the rendezvous shape. By adjusting the rendezvous time and keeping other parameters unchanged, we obtain the results listed in Table 3, where N/A means no solution. This shows that whether the rendezvous time is long or short, the proposed method always has the minimum Δ V and the minimum T a , m a x . Zeng’s method has no solution when the rendezvous time is short. This is because Zeng’s method does not consider the monotonicity of ϕ 1 and ϕ 2 but uses safety as a reference quantity to test the quality of the trajectories. Thus, there may be no solution when the rendezvous time is too short. Moreover, the safety of trajectories is considered by the proposed method in the design of the elevation angle, but it is not considered in the design of the radius. Hence, the trajectories may intersect with each other when the rendezvous time is shortened.
Let the rendezvous time be 9 years (56.5089 TU), N = 3, and the orbit elements except the orbital inclination i 2 are the same as the above example. Table 4 lists the results. This shows that the proposed method is superior to the existing two methods when the elevation angle is large. In addition, Novak’s method often fails when both the initial and target orbital inclination, i 1 and i 2 , are larger than 50 . Moreover, when i 1 = i 2 50 , the rendezvous trajectories designed by Novak’s method are not in the same plane because the projection problem of the true anomaly and azimuthal angle is not considered. Although Zeng’s method is feasible at all elevation-angles, Δ V and T a , m a x cannot be reduced.

4.2. Rendezvous Mission B from Outer Orbit to Inner Orbit

For mission B, the semi-major axes, eccentricities, and orbital inclinations of the initial and target orbit elements are set to be the same as those for mission A, respectively, which are also listed in Table 1. Let N = 3, t f = 9 years (56.5089 TU), n 3 = −20 and n 4 = −30. The trajectories are shown in Figure 6, and the results are listed in Table 5. Obviously, in this rendezvous problem, Novak’s method has no solution. The trajectory designed by Zeng’s method intersects itself, and both the velocity increment and maximal thrust acceleration are large. The proposed method is to transform the elevation angle as far as possible at the place with the largest radius, so it can be seen from Figure 7 that there will be a relatively large thrust acceleration level in the initial stage of rendezvous, but compared with Zeng’s method, its T a , m a x and Δ V are small.

4.3. Adjoint Estimation with Shape-Based Method

The velocity increment and maximal thrust acceleration of the rendezvous trajectories designed by Zeng’s method are large and intersect themselves in many cases, so they are not suitable to be used as the nominal trajectories. The following only compares the performance of the nominal trajectories designed by the proposed method with Novak’s method. Both MaxFunEvals and MaxIter options of the fsolve function in Matlab R2018a are set to 2000.
We consider mission A and make the orbit elements the same as those of mission A in Table 1 except i 2 . Let i 2 be 35 , N be 3 and t f be 7 years (43.9514TU). The spacecraft performance parameters are listed in Table 6. Here, we use the percentage of converged (POC) cases out of 100 samples to evaluate the nominal trajectories. The only difference between these 100 examples is that the initial adjoint of mass is randomly selected between 0 and the estimated value λ m t 0 [21]. Table 7 reports the POC and the initial values of the adjoint variables λ x t 0 and λ m t 0 estimated by each method with different nominal trajectories. It also shows the actual initial adjoint variables of the energy-optimal problem. The initial adjoint variables of the fuel-optimal problem are obtained by the homotopic approach based on the energy-optimal problem. Figure 8 shows the shape and normalized thrust u = T / T m a x generated by each method.
Table 7 shows that the POC of the proposed method is higher than Novak’s method, because the trajectory designed by the proposed method is closer to the energy-optimal trajectory in shape than that designed by Novak’s method; it also consumes less Δ V and has smaller T a , m a x . At the same time, compared with the non-normalized adjoint variables, the POC of the normalized adjoint variables has been significantly improved.
Table 8 presents the simulation results for cases where only the i 2 of the target orbit is changed, while other parameters are kept unchanged. It can be seen that when the adjoint variables are not normalized, the POC of Novak’s method is higher than that of the proposed method in the case of a small elevation angle, but Novak’s method is not applicable in the case of a large elevation angle. When the adjoint variables are normalized, the POC of the two methods is significantly improved, and the POC of the proposed method is significantly higher than that of Novak’s method.

5. Conclusions

This paper proposed a new shape-based method to solve 3D rendezvous problems, keeping the advantages of Novak’s method and Zeng’s method. It also improves the out-of-plane shaping functions, which makes the proposed rendezvous trajectories more flexible and more applicable. The proposed method can generate an arbitrary 3D rendezvous trajectory and can significantly reduce the velocity increment and maximal thrust acceleration for most test cases. Compared with other shape-based methods, the proposed method does not need to optimize free parameters. At the same time, this paper adds the adjoint normalization into the process of combining the shape-based method with the optimization problem and takes the proposed trajectories as the nominal orbits to provide the initial estimates of the adjoint variables for the indirect method. It is proved that the shape-based method proposed in this paper can be effectively applied to the optimization problems of orbital rendezvous problems when the inclinations of the initial and target orbits differ greatly.

Author Contributions

T.Z. and D.W. completed preliminary research and provided the numerical part; T.Z. and F.J. conceived and wrote the paper; F.J. and H.Z. supervised the overall work and reviewed the paper. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (No. 12022214).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

This Appendix shows the expressions of r , r , r , φ , φ , and φ .
The derivative and double derivative of φ with respect to θ can be expressed as
φ = e r , z e r , x 2 + e r , y 2 e r , z e r , x e r , x + e r , y e r , y e r , x 2 + e r , y 2
φ = e r , z e r , x 2 + e r , y 2 e r , z e r , x 2 + e r , x e r , x + e r , y 2 + e r , y e r , y e r , x 2 + e r , y 2 + e r , z e r , x e r , x + e r , y e r , y 2 e r , x 2 + e r , y 2 3 2
ϕ = e r , z e r , x 2 + e r , y 2 + e r , z e r , x e r , x + e r , y e r , y e r , x 2 + e r , y 2 3 e r , z e r , x e r , x + e r , y e r , y 3 e r , x 2 + e r , y 2 5 2 e r , z e r , x 2 + e r , x e r , x + e r , y 2 + e r , y e r , y + e r , z 3 e r , x e r , x + e r , x e r , x + 3 e r , y e r , y + e r , y e r , y e r , x 2 + e r , y 2 + e r , z e r , x 2 + e r , x e r , x + e r , y 2 + e r , y e r , y e r , x e r , x + e r , y e r , y + 2 + e r , z e r , x e r , x + e r , y e r , y 2 e r , x 2 + e r , y 2 3 2
where
e r , x ( y , o r , z ) = ϕ cos f 1 p 1 , x ( y , o r , z ) + sin f 1 q 1 , x ( y , o r , z ) ϕ cos f 2 p 2 , x ( y , o r , z ) + sin f 2 q 2 , x ( y , o r , z ) + ϕ p 1 , x ( y , o r , z ) f 1 sin f 1 + q 1 , x ( y , o r , z ) f 1 cos f 1 + ( 1 ϕ ) p 2 , x ( y , o r , z ) f 2 sin f 2 + q 2 , x ( y , o r , z ) f 2 cos f 2
e r , x ( y , o r , z ) = ψ cos f 1 p 1 , x ( y , o r , z ) + sin f 1 q 1 , x ( y , o r , z ) + ψ p 1 , x ( y , o r , z ) f 1 sin f 1 + q 1 , x ( y , o r , z ) f 1 cos f 1 ψ cos f 2 p 2 , x ( y , o r , z ) + sin f 2 q 2 , x ( y , o r , z ) ψ f 2 p 2 , x ( y , o r , z ) sin f 2 + f 2 q 2 , x ( y , o r , z ) cos f 2 + ψ p 1 , x ( y , o r , z ) f 1 sin f 1 + q 1 , x ( y , o r , z ) f 1 cos f 1 + ψ p 1 , x ( y , o r , z ) f 1 2 cos f 1 p 1 , x ( y , o r , z ) f 1 sin f 1 + ψ q 1 , x ( y , o r , z ) f 1 2 sin f 1 + q 1 , x ( y , o r , z ) f 1 cos f 1 + ψ p 2 , x ( y , o r , z ) f 2 sin f 2 + q 2 , x ( y , o r , z ) f 2 cos f 2 + ( 1 ψ ) p 2 , x ( y , o r , z ) f 2 2 cos f 2 p 2 , x ( y , o r , z ) f 2 sin f 2 + ( 1 ψ ) q 2 , x ( y , o r , z ) f 2 2 sin f 2 + q 2 , x ( y , o r , z ) f 2 cos f 2
e r , x ( y , o r , z ) = ϕ cos f 1 p 1 , x ( y , o r , z ) + sin f 1 q 1 , x ( y , o r , z ) ϕ cos f 2 p 2 , x ( y , o r , z ) + sin f 2 q 2 , x ( y , o r , z ) + 3 ϕ p 1 , x ( y , o r , z ) f 1 sin f 1 + q 1 , x ( y , o r , z ) f 1 cos f 1 ϕ p 2 , x ( y , o r , z ) f 2 sin f 2 + q 2 , x ( y , o r , z ) f 2 cos f 2 + 3 ϕ p 1 , x ( y , o r , z ) f 1 sin f 1 p 1 , x ( y , o r , z ) f 1 2 cos f 1 + q 1 , x ( y , o r , z ) f 1 cos f 1 q 1 , x ( y , o r , z ) f 1 2 sin f 1 ϕ f 2 p 2 , x ( y , o r , z ) sin f 2 f 2 2 p 2 , x ( y , o r , z ) cos f 2 + f 2 q 2 , x ( y , o r , z ) cos f 2 f 2 2 q 2 , x ( y , o r , z ) sin f 2 + ϕ 2 p 1 , x ( y , o r , z ) f 1 f 1 cos f 1 + p 1 , x ( y , o r , z ) f 1 3 sin f 1 p 1 , x ( y , o r , z ) f 1 sin f 1 + ϕ p 1 , x ( y , o r , z ) f 1 f 1 cos f 1 2 q 1 , x ( y , o r , z ) f 1 f 1 sin f 1 + ϕ q 1 , x ( y , o r , z ) f 1 3 cos f 1 + q 1 , x ( y , o r , z ) f 1 cos f 1 q 1 , x ( y , o r , z ) f 1 f 1 sin f 1 + ( 1 ϕ ) 2 p 2 , x ( y , o r , z ) f 2 f 2 cos f 2 + p 2 , x ( y , o r , z ) f 2 3 sin f 2 p 2 , x ( y , o r , z ) f 2 sin f 2 + ( 1 ϕ ) p 2 , x ( y , o r , z ) f 2 f 2 cos f 2 2 q 2 , x ( y , o r , z ) f 2 f 2 sin f 2 + ( 1 ϕ ) q 2 , x ( y , o r , z ) f 2 3 cos f 2 + q 2 , x ( y , o r , z ) f 2 cos f 2 q 2 , x ( y , o r , z ) f 2 f 2 sin f 2
f j = p j , x q j , y p j , y q j , x d e r f j , ( 0 ) , j = 1 , 2
f j = p j , x q j , y p j , y q j , x d e r f j , ( 1 ) d e r f j , ( 0 ) 2 , j = 1 , 2
f j = p j , x q j , y p j , y q j , x d e r f j , ( 2 ) d e r j j , ( 0 ) 2 d e r f j , ( 1 ) 2 d e r f j , ( 0 ) 3 , j = 1 , 2
d e r f j , ( 0 ) , d e r f j , ( 1 ) , d e r f j , ( 2 ) , j=1,2 are defined as
d e r f j , ( 0 ) = p j , x 2 + q j , x 2 sin θ 2 cos θ p j , x p j , y + q j , x q j , y sin θ + p j , y 2 + q j , y 2 cos 2 θ
d e r f j , ( 1 ) = p j , x 2 p j , y 2 + q j , x 2 q j , y 2 sin 2 θ 2 p j , x p j , y + q j , x q j , y cos 2 θ
d e r f j , ( 2 ) = 4 p j , x p j , y + q j , x q j , y sin 2 θ + 2 p j , x 2 p j , y 2 + q j , x 2 q j , y 2 cos 2 θ
The derivative and double derivative of the interpolation function are
if r 1 < r 2
ϕ = b + c n 1 β n 1 1 + d n 2 β n 2 1 θ 2 θ 1
ϕ = c n 1 n 1 1 β n 1 2 + d n 2 n 2 1 β n 2 2 θ 2 θ 1 2
ϕ = c n 1 n 1 1 n 1 2 β n 1 3 + d n 2 n 2 1 n 2 2 β n 2 3 θ 2 θ 1 3
if r 1 > r 2
ϕ = b + c n 3 ( β + 1 ) n 3 1 + d n 4 ( β + 1 ) n 4 1 θ 2 θ 1
ϕ = c n 3 n 3 1 ( β + 1 ) n 3 2 + d n 4 n 4 1 ( β + 1 ) n 4 2 θ 2 θ 1 2
ϕ = c n 3 n 3 1 n 3 2 β n 3 3 + d n 4 n 4 1 n 4 2 β n 4 3 θ 2 θ 1 3
Let d e r ( 0 ) denote the denominator of r, and d e r ( 1 ) , d e r ( 2 ) , and d e r ( 3 ) denote the derivative, double derivative, and triple derivative of the denominator of r with respect to f, respectively:
d e r r ( 0 ) = k 0 + k 1 f + k 2 f 2 + ( k 3 + k 4 f ) cos f + ( k 5 + k 6 f ) sin f d e r r ( 1 ) = k 1 + 2 k 2 f k 3 sin f + k 4 cos f k 4 f sin f + k 5 cos f + k 6 sin f + k 6 f cos f d e r r ( 2 ) = 2 k 2 k 3 cos f 2 k 4 sin f k 4 f cos f k 5 sin f + 2 k 6 cos f k 6 f sin f d e r r ( 3 ) = k 3 sin f 3 k 4 cos f + k 4 f sin f k 5 cos f 3 k 6 sin f k 6 f cos f
Then, r , r , and r are
r = d e r r ( 1 ) d e r r ( 0 ) 2 f
r = d e r r ( 2 ) d e r r ( 0 ) + 2 d e r r ( 1 ) 2 d e r ( 0 ) 3 f 2 + d e r r ( 1 ) d e r r ( 0 ) 2 f
r = 3 d e r r ( 2 ) d e r r ( 0 ) + 6 d e r r ( 1 ) 2 d e r r ( 0 ) 3 f f + d e r r ( 1 ) d e r r ( 0 ) 2 f + d e r r ( 3 ) d e r r ( 0 ) 2 + 6 d e r r ( 0 ) d e r r ( 1 ) d e r r ( 2 ) 6 d e r r ( 1 ) 3 d e r r ( 0 ) 4 f 3
where
f = p x q y p y q x d e r f ( 0 )
f = p x q y p y q x d e r f ( 1 ) d e r f ( 0 ) 2
f = p x q y p y q x d e r f ( 2 ) d e r f ( 0 ) 2 d e r f ( 1 ) 2 d e r f ( 0 ) 3
d e r f ( 0 ) , d e r f ( 1 ) , and d e r f ( 2 ) are denoted as
d e r f ( 0 ) = p x 2 + q x 2 sin θ 2 cos θ p x p y + q x q y sin θ + p y 2 + q y 2 cos 2 θ
d e r f ( 1 ) = p x 2 p y 2 + q x 2 q y 2 sin 2 θ 2 p x p y + q x q y cos 2 θ
d e r f ( 2 ) = 4 p x p y + q x q y sin 2 θ + 2 p x 2 p y 2 + q x 2 q y 2 cos 2 θ

References

  1. Betts, J. Survey of Numerical Methods for Trajectory Optimization. J. Guid. Control. Dyn. 1998, 21, 193–207. [Google Scholar] [CrossRef]
  2. David, G.H. Conversion of Optimal Control Problems into Parameter Optimization Problems. J. Guid. Control. Dyn. 1997, 20, 57–60. [Google Scholar]
  3. Caruso, A.; Quarta, A.A.; Mengali, G. Comparison between Direct and Indirect Approach to Solar Sail Circle-to-circle Orbit Raising Optimization. Astrodynamics 2019, 3, 273–284. [Google Scholar] [CrossRef]
  4. Yan, H.; Wu, H. Initial Adjoint-Variable Guess Technique and Its Application in Optimal Orbital Transfer. J. Guid. Control. Dyn. 1999, 22, 490–492. [Google Scholar] [CrossRef]
  5. Mantia, M.; Casalino, L. Indirect Optimization of Low-Thrust Capture Trajectories. J. Guid. Control. Dyn. 2006, 29, 1011–1014. [Google Scholar] [CrossRef]
  6. Martell, C.A.; Lawton, J.A. Adjoint Variable Solutions via an Auxiliary Optimization Problem. J. Guid. Control. Dyn. 2015, 18, 1267–1272. [Google Scholar] [CrossRef]
  7. Kluever, C.A.; Pierson, B.L. Optimal Low-Thrust Earth-Moon Transfers with a Switching Function Structure. J. Astronaut. Sci. 1994, 42, 269–283. [Google Scholar]
  8. Kluever, C.A.; Pierson, B.L. Optimal Low-thrust Three-Dimensional Earth-Moon Trajectories. J. Guid. Control Dyn. 1995, 18, 830–837. [Google Scholar] [CrossRef]
  9. Gao, Y.; Kluever, C. Low-Thrust Interplanetary Orbit Transfers Using Hybrid Trajectory Optimization Method with Multiple Shooting. In Proceedings of the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Providence, RI, USA, 16–19 August 2004. [Google Scholar]
  10. Yang, H.; Tang, G.; Jiang, F. Optimization of Observing Sequence Based on Nominal Trajectories of Symmetric Observing Configuration. Astrodynamics 2018, 107, 25–37. [Google Scholar] [CrossRef]
  11. Petropoulos, A.E.; Longuski, J.M. Shape-Based Algorithm for the Automated Design of Low-Thrust, Gravity Assist Trajectories. J. Spacecr. Rocket. 2004, 41, 787–796. [Google Scholar] [CrossRef] [Green Version]
  12. Wall, B.J.; Conway, B.A. Shape-Based Approach to Low-Thrust Rendezvous Trajectory Design. J. Guid. Control. Dyn. 2009, 32, 95–101. [Google Scholar] [CrossRef]
  13. Wall, B.J.; Pols, B.; Lanktree, B. Shape-based Approximation Method for Low-thrust Interception and Rendezvous Trajectory Design. Adv. Astronaut. Sci. 2010, 136, 1447–1458. [Google Scholar]
  14. Wall, B.J. Shape-Based Approximation Method for Low-Thrust Trajectory Optimization. In Proceedings of the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, HI, USA, 18–21 August 2008. [Google Scholar]
  15. Novak, D.M.; Vasile, M. Improved Shaping Approach to the Preliminary Design of Low-Thrust Trajectories. J. Guid. Control. Dyn. 2011, 34, 128–147. [Google Scholar] [CrossRef] [Green Version]
  16. Xie, C.; Zhang, G.; Zhang, Y. Simple Shaping Approximation for Low-Thrust Trajectories Between Coplanar Elliptical Orbits. J. Guid. Control. Dyn. 2015, 38, 2448–2455. [Google Scholar] [CrossRef]
  17. Xie, C.; Zhang, G.; Zhang, Y. Shaping Approximation for Low-Thrust Trajectories with Large Out-of-Plane Motion. J. Guid. Control Dyn. 2016, 39, 2776–2786. [Google Scholar] [CrossRef]
  18. Zeng, K.; Geng, Y.; Wu, B. Shape-Based Analytic Safe Trajectory Design for Spacecraft Equipped with Low-Thrust Engines. Aerosp. Sci. Technol. 2017, 62, 87–97. [Google Scholar] [CrossRef]
  19. Jiang, F.; Baoyin, H.; Li, J. Practical Techniques for Low-Thrust Trajectory Optimization with Homotopic Approach. J. Guid. Control. Dyn. 2012, 35, 245–258. [Google Scholar] [CrossRef]
  20. Petukhov, V.G. Optimization of Interplanetary Trajectories for Spacecraft with Ideally Regulated Engines Using the Continuation Method. Cosm. Res. 2008, 46, 219–232. [Google Scholar] [CrossRef]
  21. Jiang, F.; Tang, G.; Li, J. Improving Low-Thrust Trajectory Optimization by Adjoint Estimation with Shape-Based Path. J. Guid. Control. Dyn. 2017, 40, 1–8. [Google Scholar] [CrossRef]
  22. Taheri, E.; Kolmanovsky, I.; Atkins, E. Enhanced Smoothing Technique for Indirect Optimization of Minimum-Fuel Low-Thrust Trajectories. J. Guid. Control. Dyn. 2016, 39, 2500–2511. [Google Scholar] [CrossRef] [Green Version]
  23. Lewallen, J.M.; Schwausch, O.A.; Tapley, B.D. Coordinate System Influence on the Regularized Trajectory Optimization Problem. J. Spacecr. Rocket. 1971, 8, 15–20. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Illustration of the spherical coordinates.
Figure 1. Illustration of the spherical coordinates.
Aerospace 08 00315 g001
Figure 2. Diagram of the relationship between f and θ .
Figure 2. Diagram of the relationship between f and θ .
Aerospace 08 00315 g002
Figure 3. Effect of the shaping parameters on ϕ .
Figure 3. Effect of the shaping parameters on ϕ .
Aerospace 08 00315 g003
Figure 4. Rendezvous results of three methods when r 1 < r 2 .
Figure 4. Rendezvous results of three methods when r 1 < r 2 .
Aerospace 08 00315 g004
Figure 5. Thrust acceleration T a of three methods when r 1 < r 2 .
Figure 5. Thrust acceleration T a of three methods when r 1 < r 2 .
Aerospace 08 00315 g005
Figure 6. Rendezvous result of three methods when r 1 > r 2 .
Figure 6. Rendezvous result of three methods when r 1 > r 2 .
Aerospace 08 00315 g006
Figure 7. Thrust acceleration T a of each method when r 1 > r 2 .
Figure 7. Thrust acceleration T a of each method when r 1 > r 2 .
Aerospace 08 00315 g007
Figure 8. Rendezvous result of three methods.
Figure 8. Rendezvous result of three methods.
Aerospace 08 00315 g008
Table 1. Orbit elements of two missions.
Table 1. Orbit elements of two missions.
Mission A ( r 1 < r 2 ) Mission B ( r 1 > r 2 )
Orbit ElementsInitial OrbitTarget OrbitInitial OrbitTarget Orbit
Semi-major axis1441
Eccentricity0.010.10.10.01
Inclination 5 65 65 5
Right ascension 10 10 10 10
Argument of periapsis 10 10 10 10
True anomaly 0 100 0 100
Table 2. Orbit rendezvous results when r 1 < r 2 .
Table 2. Orbit rendezvous results when r 1 < r 2 .
Method Δ V (VU) T a , max ( LU / TU 2 )
Proposed1.59380.1527
Novak2.80700.2943
Zeng10.24393.4581
Table 3. Orbit rendezvous results with different t f when r 1 < r 2 .
Table 3. Orbit rendezvous results with different t f when r 1 < r 2 .
Proposed Novak Zeng
t f Δ V (VU) T a , max ( LU / TU 2 ) Δ V (VU) T a , max ( LU / TU 2 ) Δ V (VU) T a , max ( LU / TU 2 )
901.77640.2229N/AN/A2.44250.2997
801.70600.2040N/AN/A2.37230.2442
701.64740.1831N/AN/A2.72640.3755
601.60480.16072.73950.27976.42891.2344
501.58380.13783.00000.5310N/AN/A
401.72490.2099N/AN/AN/AN/A
302.34490.4569N/AN/AN/AN/A
Table 4. Orbit rendezvous results with different i 2 when r 1 < r 2 .
Table 4. Orbit rendezvous results with different i 2 when r 1 < r 2 .
Proposed Novak Zeng
i 2 Δ V (VU) T a , max ( LU / TU 2 ) Δ V (VU) T a , max ( LU / TU 2 ) Δ V (VU) T a , max ( LU / TU 2 )
75 2.67850.7977N/AN/A9.49522.4343
65 1.59380.15272.80700.294310.24393.4581
55 1.19060.06511.88300.100610.25244.1412
45 0.95530.03811.28620.07069.85554.6166
35 0.78870.02570.89310.05489.25355.0227
25 0.65390.02280.65820.03938.62345.7724
15 0.54650.02120.53840.02568.13947.2918
Table 5. Orbit rendezvous results when r 1 > r 2 .
Table 5. Orbit rendezvous results when r 1 > r 2 .
Method Δ V (VU) T a , max ( LU / TU 2 )
Proposed1.83450.3646
NovakN/AN/A
Zeng14.25086.4255
Table 6. Spacecraft performance parameters.
Table 6. Spacecraft performance parameters.
I sp ( s ) T max ( N ) m 0 ( kg )
30000.85000
Table 7. Initial costates and POC of different methods.
Table 7. Initial costates and POC of different methods.
Method λ x λ m POC
Proposed 0.8032 , 0.0005 , 0.1797 , 0.1241 , 0.6611 , 0.8892 T 0.8963 5 %
Proposed (normalization) 0.4162 , 0.0003 , 0.0931 , 0.0643 , 0.3426 , 0.4608 T 0.4645 69 %
Novak 0.9528 , 0.0155 , 0.1093 , 0.0433 , 1.1361 , 0.7288 T 1.2450 3 %
Novak (normalization) 0.4141 , 0.0067 , 0.0475 , 0.0188 , 0.4938 , 0.3168 T 0.5411 23 %
Energy-optimal 0.2012 , 0.0012 , 0.1084 , 0.0020 0.1083 , 0.2867 T 0.4779 N/A
Fuel-optimal 0.2506 , 0.0000 , 0.0841 , 0.0238 0.1397 , 0.2158 T 0.5017 N/A
Table 8. POC of different methods with different i 2 .
Table 8. POC of different methods with different i 2 .
i 2 ProposedProposed (Normalization)NovakNovak (Normalization)
75 0 % 33 % N/AN/A
65 0 % 47 % N/AN/A
55 0 % 11 % 0 % 0 %
45 0 % 42 % 0 % 24 %
35 5 % 69 % 3 % 23 %
25 19 % 87 % 30 % 51 %
15 68 % 99 % 99 % 100 %
5 100 % 100 % 96 % 100 %
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, T.; Wu, D.; Jiang, F.; Zhou, H. A New 3D Shaping Method for Low-Thrust Trajectories between Non-Intersect Orbits. Aerospace 2021, 8, 315. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace8110315

AMA Style

Zhang T, Wu D, Jiang F, Zhou H. A New 3D Shaping Method for Low-Thrust Trajectories between Non-Intersect Orbits. Aerospace. 2021; 8(11):315. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace8110315

Chicago/Turabian Style

Zhang, Tongxin, Di Wu, Fanghua Jiang, and Hong Zhou. 2021. "A New 3D Shaping Method for Low-Thrust Trajectories between Non-Intersect Orbits" Aerospace 8, no. 11: 315. https://0-doi-org.brum.beds.ac.uk/10.3390/aerospace8110315

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop