International Journal of Aerospace Engineering

Volume 2019, Article ID 1364834, 18 pages

https://doi.org/10.1155/2019/1364834

## Low-Thrust Trajectory Design Using Finite Fourier Series Approximation of Pseudoequinoctial Elements

National University of Defense Technology, Changsha, Hunan 410073, China

Correspondence should be addressed to Haiyang Li; nc.ude.tdun@gnayiahil

Received 17 March 2019; Revised 25 June 2019; Accepted 3 September 2019; Published 4 November 2019

Academic Editor: Seid H. Pourtakdoust

Copyright © 2019 Wanmeng Zhou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A new low-thrust trajectory design method is proposed that is based on the finite Fourier series method with pseudoequinoctial elements rather than the more common cylindrical coordinate components. The bijection relation between the elements and control variables is ensured by introducing an additional equality constraint derived from the angular momentum conservation. The guidance law and on-line control variables are obtained by applying inverse dynamics and the framework of inverse simulation technology, respectively. The pseudoequinoctial finite Fourier series method has the advantages of both the Fourier series and the perturbation analysis methods. For two-body problems, three cases were studied: the Earth to Mars, 1989ML, and Tempel-1 missions. Regarding the design of a rendezvous trajectory with a large inclination angle and a high eccentricity rate, this method yields a broader range of feasible results than the traditional Fourier series method. The circular restricted three-body problem was solved for the first time using the pseudoequinoctial finite Fourier series method combined with the patched conics method. The low-thrust Earth-Moon transfer was analyzed, and the results show that this method improves window analysis efficiency and guarantees precision of the initial geocentric trajectory for the low-thrust transfer.

#### 1. Introduction

In recent years, low-thrust spacecraft have increasingly attracted interest in the field of engineering because of their long-lasting and fuel-saving features. These include the National Aeronautics and Space Administration’s (NASA) Deep Space 1 [1] and Dawn [2], the European Space Agency’s (ESA) SMART-1 [3], and BepiColombo [4], a joint mission between ESA and Japan Aerospace Exploration Agency (JAXA). The trajectory design of these low-thrust vehicles has also attracted much attention. In a global trajectory optimization competition (GTOC 2), the trajectory design problem for low-thrust vehicles was proposed to explore asteroids and issues associated with them [5]. The challenges in the search for the optimal trajectory of a low-thrust, deep-space exploration vehicle are greater than those of an impulsive-thrust design due to features such as the long time of flight (TOF) and multiple solutions of the feasible flight trajectory. Given the initial and final position, velocity, and TOF, the thrust profile can be determined by optimal fuel or energy under maximum thrust amplitude, which is also called a two-point boundary value (TPBV) problem [6, 7].

There are direct and indirect methods for solving TPBV problems. In indirect methods [8], according to Pontryagin’s maximum principle, the derivative of the Hamiltonian function with respect to the control input equals zero. The control input can then be represented by both state and costate variables [9]. Substituting the control back into the differential equations of the costate and state variables, the simultaneous equations are then solved based on the boundary and transversal conditions. The disadvantage of these methods is that it is difficult to assign the costates’ appropriate initial values when solving the problem, for their unapparent physical significances. In direct methods, such as point collocation [10] and the pseudospectral method [11], either the control variables or both the state and control variables are discretized to form a set of nonlinear equations with unknown parameters that are solved by nonlinear programming algorithms. Although these methods do not require the introduction of costate variables, the discrete points may become increasingly dense during the iteration process, which makes the calculation time consuming and potentially unstable. Other variants of direct methods are the shape-based methods, in which the flight states are suitably represented by parameterized nonlinear functions with certain assumptions and use inverse dynamics to retrieve the control corresponding to the assumed functional form. These methods require fewer collocation points and design parameters to obtain higher computational efficiency. Based on previous discussions of the optimal direction of a fixed thrust [12–14], early shape-based methods on two-dimensional problems were generally based on the tangential thrust assumption and the analytic results were typically expressed in polar coordinates [15–17] such as in the logarithmic spiral method, Forbes spiral method, exponential sinusoid method, and inverse polynomial method (IP). As the number of shape parameters increases, the number of constraint conditions that can be satisfied also increases. The design parameters can be derived directly from the constraint conditions without numerical iterations. To solve three-dimensional (3D) problems, there are several methods such as IP in cylindrical coordinates [18], the pseudoequinoctial method (PE) [19] and the spherical shaping method [20, 21] on the basis of the perturbation functions, and the hodograph method which considers the velocity variable as the parameterized shape [22]. The IP method cannot ensure the proper amplitude of the thrust, while PE and spherical shaping methods require additional design parameters to satisfy the TOF or thrust constraints [20]; neither can these methods provide the motion states directly related to the transfer time. For the hodograph method, an appropriate combination of basic functions should be selected first to obtain the optimal velocity.

The finite Fourier series method (FFS) is a new method to overcome these deficiencies. It applies time as an independent variable and represents the position variables with an orthogonal polynomial with a finite number of trigonometric terms, i.e., Fourier series [23–25]. The orthogonal polynomials have been proven to solve optimal control problems successfully [26, 27]. The Fourier series in particular can represent any periodical function and could even theoretically represent the switching-type function with sufficient trigonometric function terms [28], which introduces more shape categories into the trajectory design. On the other hand, most research on the shape-based methods is applied on two-body problems, but few are applied for circular restricted three-body problems (CRTBPs). FFS combined with the simplified analysis method is the only shape-based method proposed so far to solve CRTBPs [29].

Although FFS produces excellent results in transfer trajectory design between coplanar circular orbits [30], the increasing eccentricity and inclination angle of the transfer orbit will cause the total velocity increment to gradually diverge from the near-optimum solution and there will be more infeasible points in the launch window graph [31]. The first major reason is because there are a fixed number of terms selected which leads to a relatively fixed fitting shape. Another reason is because the expression of 3D-FFS in a cylindrical coordinate system has limited approximation capability. During the long-time transfer, the motion in the -direction fluctuates periodically due to changing eccentricity, inclination angle, and the semimajor axis, which is difficult to match. In this study, the PE elements were introduced into the modified scheme of the FFS, denoted as PE-FFS. The PEs are similar to modified equinoctial elements (MEEs), which can provide smooth variations of instantaneous motion states compared to the -axis components of a cylindrical coordinate system [32]. For two-body problems, different numbers of terms in each PE are first chosen for the different relative positions of the initial and final radius vectors. A feasible shape library expressed by combinations of these terms is then used to produce a wider range of feasible solutions. For CRTBPs, the new method proposed is the PE-FFS combined with the patched conics method. The capture situation and feasibility of the transfer trajectory are comprehensively analyzed through the state change at the capture point of the lunar sphere of influence (LSOI).

This paper is organized as follows: First, the dynamic equations in the cylindrical coordinate system and the Gaussian perturbation equations of the modified equinoctial orbit elements are established separately under the two-body assumption. Then, the inverse dynamic model and inverse simulation technique are proposed to give the relationship between the state and control variables, and an additional condition is proposed to form the PE-FFS method. The third section presents the results of studies of the two-body rendezvous cases, which include transfers to Mars, an asteroid, and a comet. The results from the PE-FFS method are compared with those from the other methods. Finally, the CRTBP for the Earth-Moon transfer is investigated. The quick trajectory design and the launch window analysis for the Earth-Moon mission are realized by combining the PE-FFS and patched conics method.

#### 2. Pseudoequinoctial Finite Fourier Series Method

##### 2.1. Dynamic Modeling of a Two-Body System

A typical two-body problem [19] can be described as where refers to the gravitational constant of the central body and refers to the acceleration produced by the external force. In traditional shape-based methods, the state variable is represented by the polar radius , the phase angle , and the height along the -axis. The components of equation (1) [32, 33] can be expressed as where and refer to the acceleration components in the local cylindrical coordinates. The definition of directions in the cylindrical coordinate system is shown in Figure 1.