Abstract
Space trajectory design is usually addressed as an optimal control problem. Although it relies on the classic theory of optimal control, this branch possesses some peculiarities that led to the development of ad hoc techniques, which can be grouped into two categories: direct and indirect methods. This paper gives an overview of the principal techniques belonging to the direct methods. The technique known as “direct transcription and collocation” is illustrated by considering Hermite-Simpson, high-order Gauss-Lobatto, and pseudospectral methods. Practical examples are given, and several hints to improve efficiency and robustness are implemented.
1. Introduction
In January 1959, the Soviet Union launched Luna 1, a lunar probe that became the first spacecraft placed in heliocentric orbit (http://en.wikipedia.org/wiki/Luna_1). This event marked a new era for deep space exploration. In over half a century, various probes visited the major planets of the Solar System, and they provided a remarkable scientific return. Nowadays, many space agencies have established ambitious space programs, whose accomplishment stimulated the development of new technologies.
Electric propulsion is one of those technologies that improve the efficiency of space transport. Electric motors have a specific impulse that is approximately ten times that of the chemical engines, and therefore they allow a considerable saving of propellant mass, which in turn makes it possible to embark heavier instruments, so increasing the scientific return of a mission. Electric propulsion also involves long working hours, extended launch windows, and flexible and precise control abilities. Electric propulsion has been used in NASA’s Deep Space 1 in 1998 [1], in JAXA’s Hayabusa in 2003 [2], and in ESA’s SMART-1 in 2003 [3]. Many future space missions foresee the use of this technology. If compared to that of the chemical engines, electric propulsion produces low levels of thrust, which cannot be modelled as producing instantaneous velocity changes [4]. This is also the case of solar sails, an emerging and prominent technology demonstrated by the recent success of JAXA’s IKAROS [5] and NASA’s NanoSail-D2 [6]. Solar sails are not propellant-constrained, and the effect of the solar radiation pressure on a large surface can be modelled as producing a low-thrust acceleration [7, 8]. The combination of low-thrust propulsion with gravity assists is among the most promising techniques for deep space explorations, which can be utilized to save fuel and shorten the flight time [9–12].
Low-thrust space trajectories are studied as a specialization of the optimal control problem for continuous time systems. Unfortunately, no analytic solutions exist for this problem even under the simple two-body dynamics. Thus, numerical methods must be used. Low-thrust trajectory optimization involves determining the control law (thrust magnitude and direction) and the associated transfer orbit while minimizing a given performance index (propellant mass or time-of-flight) and satisfying boundary conditions (departure and arrival orbits), midpoint conditions (at patching points), and path constraints (thrust saturation). In this branch of optimal control, two main solution methods have been devised. Within direct methods, the state and control variables are discretized, and the optimal control problem is converted into a nonlinear programming (NLP) problem [13, 14]. This process is called direct transcription. Direct methods are generally robust and can easily accommodate path constraints, but they often require much computational effort especially for multispiral trajectories. Indirect methods rely on the calculus of variations. The necessary conditions of optimality require the solution of a two-point boundary value problem (TPBVP). This method ensures rapid convergence of good starting guesses, but most of the difficulties are related to the small convergence radius, the high sensitivity to the initial costates, and their lack of physical meaning.
This paper complements previous surveys on optimal control of space trajectories [15–18] by focussing on direct methods. In particular, the direct transcription and collocation is studied. The collocation schemes considered implementing different integration methods (Hermite-Simpson, high-order Gauss-Lobatto, and pseudospectral methods). These techniques are implemented to solve practical examples, and the behavior of each method is assessed against the implementation of several issues aimed at improving the computational efficiency and the methods robustness.
The remainder of the paper is organized as follows. In Section 2, the optimal control problem in space flight mechanics is stated, and the NLP problem as well as the direct transcription concept is recalled. In Section 3 direct transcription and collocation is treated. Computational issues are discussed in Section 4 and practical cases are solved in Section 5. Concluding remarks are given in Section 6.
2. The Optimal Control Problem in Space Flight Mechanics
Given a set of first-order differential equations where is a vector of state variables, denotes vector of control variables, and represents the independent time variable, , the following performance index has to be minimized while satisfying -dimensional final boundary conditions The solution to this problem is derived by the calculus of variations, which leads to the derivation of the Euler-Lagrange equations. Let be a -dimensional constant vector of multipliers of the final boundary constraints, and let be the -dimensional variable vector of adjoint or costate multipliers of the dynamics. The augmented performance index is defined as The augmented performance index (4) embeds the dynamics (1) as well as the final boundary conditions (3). The problem consists in deriving the necessary conditions for a stationary point of [19]. This is achieved by imposing that its first variation is zero, namely, . In order to write the necessary conditions in a compact form, it is convenient to define the Hamiltonian The necessary conditions for optimality, also referred to as the Euler-Lagrange equations [20, 21], are where the subscript denotes partial derivation. The first of (6) is equivalent to (1), the second describes the dynamics of the costates, and the third is an algebraic equation for the control functions. This differential-algebraic system must be solved together with the final boundary conditions (3) and the following transversality conditions: With given, the problem represents a two-point boundary value problem. The last of (6) is an application of the Pontryagin maximum principle [22]. A more general expression is in fact where defines the domain of feasible controls. The maximum principle states that the control variables must be chosen to optimize the Hamiltonian at every instant of time: the solution of the optimal control problem is an extremum for . In essence, the maximum principle is a constrained optimization problem in the function at all values of .
2.1. The Optimal Trajectory Design Problem
An optimal trajectory design problem is a specialization of the classical optimal control problem above. In space flight mechanics, the equations of motion have the following form [16]: In (9), and are the spacecraft position and velocity vectors, respectively, is the gravitational vector field, is the thrust acceleration magnitude, and is the thrust direction unit vector. The control variables are and . The control acceleration magnitude is upper bounded for technological reasons; that is, . To minimize the total velocity change, and therefore the propellant mass, the objective function (2) is such that , . The Hamiltonian (5) is, therefore, where and are the costate vectors associated to the position and velocity vectors, respectively. As is bounded, the last of (6) cannot be applied, and (8) has to be used instead. Since the Hamiltonian has to be minimized at any time, the following observations can be made.(1)The thrust unit vector has to be parallel and opposite to . That is, with ; this explains why is also referred to as the primer vector.(2)The thrust magnitude has to be chosen according to the sign of the switching function: In particular, when and when . This makes the function have a bang-bang structure; that is, it is piece-wise discontinuous and it is either zero or maximum. This property is of great importance to evaluate the optimal control profile a posteriori.
The necessary conditions only guarantee that the optimal trajectory is an extremum for the Hamiltonian. Thus, to assess the optimality of the solutions, one is supposed to check the second-order conditions. However, due to the nature of the space trajectory problem, there is no upper bound to the propellant that can be consumed in one trajectory, so one may be confident that a solution that satisfies the necessary conditions is a local minimum and not a local maximum [16].
In trajectory optimization problems, the initial state, , is generally given. If the initial costate, , was given as well, the optimal solution could be obtained by integrating the Euler-Lagrange equations with an implicit “bang-bang” thrusting structure [23]. Thus, the low-thrust optimal trajectory design can be converted into a TPBVP, which consists in finding the unknown initial costate vector. This is the essence of indirect methods. Another philosophy consists instead in translating the continuous optimal control problem into a NLP problem and solving for a finite set of variables [15, 24–27]. This procedure is the direct transcription and the approach is said direct method.
2.1.1. The General Optimal Trajectory Design Problem
It is convenient to state the optimal trajectory design problem in a more general fashion, which copes with the direct approach. The dynamics are let to incorporate a number of constant parameters ; that is, and initial and final conditions can be defined within some prescribed lower and upper bounds In addition, the solution can be subject to path constraints of the form as well as simple bounds on the state variables and on the control variables The basic problem is to determine the control vectors to minimize the performance index which is written in the Mayer form [15].
2.2. The Nonlinear Programming Problem
Essentially, any numerical method for solving the trajectory optimization problem incorporates some type of Newton method to solve for a finite set of unknowns. In Section 3 it is shown how an optimal control problem can be transformed into a NLP problem [16, 27]. A NLP problem is a decisional problem concerning a scalar objective function and a vector of constraints. As opposite to the optimal control problem, no dynamics is involved in a NLP problem. Suppose that the variables must be chosen to solve subject to the equality constraints where . The Lagrangian of this problem is which is a scalar function of the variables and the Lagrange multipliers . The necessary conditions for a point () to be a constrained optimum require solving the following system: where and are the gradient of the objective function and the Jacobian of the equality constraint vector, respectively. The system (21) can be solved via a Newton method to find the variables . Given a generic initial guess , its corrections to construct the new solution () are given by solving the linear system also referred to as Karush-Kuhn-Tucker (KKT) system. In (22), is the Hessian of (20) in ; namely, It is important to observe that an equivalent way to define the search direction is to minimize the quadratic form subject to the linear constraints This is the reason why this problem is also referred to as a quadratic programming (QP) problem.
The NLP problem formulated above can be generalized to the case that occur when inequality constraints are imposed; the constraints are of the form Constraints that are strictly satisfied, that is, those for which , are called inactive; the remaining active set of constraints are on their bounds; that is, . If the active set of constraints is known, the inactive constraints are ignored and the problem is simply solved using the method for an equality constrained problem discussed above.
In summary, the general NLP problem requires finding the vectors to solve subject to the constraints and bounds In this formulation equality constraints can be imposed by setting .
2.3. Direct Transcription
With a direct approach, the solution to the optimal control problem is strictly connected to the numerical integration of the differential equations. The core of this method consists in the way the dynamics are handled; the set of differential equations governing the motion of a spacecraft can be transcribed into a finite set of equality constraints. If these are respected, then the original optimal trajectory design problem is solved within the degree of accuracy of the numerical scheme used [16, 18, 27–30].
In problem (12)–(17), the time domain can be uniformly discretized as where the time labels are referred to as mesh points or nodes; and are the initial and final time, respectively, and is the fixed step size of the discretization (a nonuniform time discretization is also possible). The states and the controls can be discretized over the mesh (30) by defining and . The discretized states and the controls are now ready to be treated as a set of NLP variables. The whole variable vector of the problem is The differential equations are replaced by a finite set of defects constraints derived by the numerical integration scheme. For instance, if a forward Euler scheme is used, the defects are of the form where . As a result of the transcription, the optimal control constraints (13)-(14) are replaced by the NLP constraints where with , , and
The first equality constraints in (34) require that the defect vectors , , are zero, thereby they satisfy the differential equations (12) within the accuracy of the numerical integration. The boundary conditions (13) are enforced directly by the equality constraints on and , and the nonlinear path constraints (14) are imposed at the grid points. In a similar fashion the objective function, either in the form (2) or (17), can be written in terms of ; namely, . The optimal trajectory design problem is so translated into the forms (27)–(29) and it can be solved as a standard NLP problem through (18)–(23).
In this derivation, the simple forward Euler scheme is used. More accurate methods are likely to be used for practical applications. This is done to keep within reasonable values, avoiding out-of-memory problems. The integration schemes influence the robustness of the method and the solution accuracy. In Section 3, several implicit quadrature methods are shown.
3. Direct Transcription and Collocation
Collocation is used to transcribe differential dynamic constraints into a set of algebraic constraints. The basic idea is to choose a polynomial up to a certain degree with a number of points in the time domain (collocation points), and to enforce the polynomials to satisfy the equations of motion at the collocation points. The peculiarity of each collocation method relies on the way the state and control variables are discretized and how the dynamic constraints are satisfied [17].
3.1. Hermite-Simpson Method
A basic form of collocation is the Hermite-Simpson method [24], illustrated in Figure 1. For each of the segments , the two end points, denoted as “nodes” (blue dots), represent the corresponding state and control NLP variables; that is, (scalar states and controls are used from now on to ease the notation). The dynamics are used to provide time derivative values at the two nodes, so the four pieces of information can be used to construct a third-order Hermite interpolate polynomial. This interpolate polynomial cannot satisfy the equations of motion at any time within , because it does that only at the nodes. Let be the state and control at , the middle point of ; this is called “collocation point” (red diamond). Enforcing makes it possible to have a polynomial that not only satisfies the dynamics at the two nodes but also does that at the collocation point. If a large number of intervals are used, the state motion approaches the real dynamics within the whole time domain.
The detailed procedure can be derived as follows [16]. Let the state function be represented on each segment through a cubic polynomial of the form which yields where are the coefficients of the polynomial. In order to simplify the argument, the time domain is transformed such that ( is the time interval of the segment). Let , , , and . Evaluating (36)-(37) at and yields which allows us to compute the four coefficients of interpolate polynomial Substituting into (36)-(37) allows us to compute the collocation point as as well as its time derivative The control variable at the collocation point can be computed by simple linear interpolation; that is, The difference between the interpolated and calculated derivatives at the collocation point defines the integration defect The NLP solver will select to drive to zero and in this way the interpolating polynomial will approximate the true dynamics within the accuracy of the numerical integration. Note that the last row of (43) is actually an implicit Hermite integration. Thus, if the collocation constraints are satisfied, the system is said to be “implicitly” integrated.
3.2. High Order Gauss-Lobatto Method
In the Hermite-Simpson method, only two nodes are used to construct a third-order polynomial. As a general rule of thumb, a lower number of segments can be handled if higher order integration is performed. Herman and Conway [28] demonstrated that higher order Gauss-Lobatto methods are more robust and more efficient than the lower-order Hermite-Simpson scheme. Figure 2 illustrates the collocation constraints of fifth-order Gauss-Lobatto method. Similar to Hermite-Simpson method, in each segment, six pieces of information of nodes are used to construct a fifth-order Hermite polynomial to approximate the state time history. Then the resulting interpolation polynomial is used to evaluate the states at the remaining two collocation points.
In the fifth-order method, the collocation points are [28] and the two collocation constraints defects to zero are
A detailed derivation of the formulas above is required when the argument is extended to arbitrary higher orders. In [31], an alternative framework for unifying arbitrary higher-order methods is presented. In this approach, the Legendre-Gauss-Lobatto (LGL) discrete points are used to improve both interpolating precision and quadrature performance. Figure 3 graphically illustrates the location of three, five, and seven LGL points. Table 1 lists the corresponding position of these LGL points (the time interval is ).
With reference to Figure 3, the blue points are the generic nodes, which are defined by while the red diamonds are the generic collocation points, denoted as such that there is one collocation point between every two adjacent nodes [31]. The nodes are used for constructing the interpolation polynomial, while the collocation points are used to formulate the defect constraints (it is noted that nodes and collocation points should not overlap in Gauss-Lobatto method). With three points, the Gauss-Lobatto method degenerates to the Hermite-Simpson method. The state in the th subinterval is approximated by the th degree Hermite interpolating polynomial where the coefficients of Hermite interpolating polynomial are determined by using the values of the states and vector field at the points ; that is, Let (49) be written as . The term depends only on the number and location of node points, whereas is a concatenated vector of state and vector field values. Thus, given a vector of states, the coefficients of the Hermite interpolation in the th subinterval, , can be found by . The values of the states located at collocation points are then The derivatives of the Hermite interpolating polynomial located at the collocation points are given by In this form the matrices and are constants, thus a system of constraints per interval is obtained as follows: and the analytic Jacobian for the defect constraints can be derived by chain rule as If the final transfer time is a variable, within each segment, the analytic Jacobian of is written as A graphical illustration of the Jacobian structure with free final time is given in Figure 4. It can be seen that each segment has the same Jacobian module; thus (53)-(54) can be calculated offline and stored before solving the NLP. In this way, the computational time can be considerably reduced compared to the use of finite difference approximation.
3.3. Pseudospectral Method
Recently, new direct methods have been applied to low-thrust trajectory optimization, the most popular one being the pseudospectral scheme [32–35]. The major difference between Gauss-Lobatto and pseudospectral collocation schemes is the way in which the interpolation polynomial is constructed and the defect constraints are defined. When the pseudospectral method is employed, a global Lagrange interpolation polynomial is constructed to approximate the state profile. The polynomial is then differentiated and evaluated at all the nodes to compute interpolated derivative values of the states, and equations of motion are used to provide physical time derivatives. The differences between the two sets of time derivatives form the defects. The formulation of defect constraints is sketched in Figure 5.
The defect constraints of the pseudospectral method can be derived as follows. The state curve is approximated by using values of at the discrete time points and the corresponding Lagrange interpolating polynomials , where are defined as Then, the derivative of is given: The time derivative of the polynomial can be expressed in compact form through where is differentiation matrix. This only depends on the chosen node spacing, so it can be calculated offline and stored. Then the defect constraints can be written as A variety of pseudospectral schemes emerged and the major differences are the way in which the nodes are selected. The most general form is Legendre-Gauss-Lobatto discretization [36], but Legendre-Gauss [34, 37] and Legendre-Gauss-Radau [35] are also commonly used. Figure 6 illustrates three different kinds of Gauss points. Table 2 lists the position of equidistant points with three kinds of Gaussian points. We can see that LG points are located on the open interval , in contrast, LGR points are located on half-open interval , and LGL points are located on closed interval .
The main characteristic of the Gauss-Lobatto method is the combination of reasonable accuracy with highly sparse constraint Jacobians and Hessians. The pseudospectral method has a more elegant form of Jacobian and it offers spectral accuracy for smooth problems, but the constraint Jacobian is much denser. So there is a balance between accuracy and efficiency. In general, the Gauss-Lobatto method use a limited number of nodes per segment () with many segments, while pseudospectral methods use many nodes per segment () with a limited number of segments.
3.4. Control Parameterization Method
The control parameterization method is a further direct transcription technique for solving optimal control problems [38–40]. In this method, only the controls are discretized and represented by a linear combination of basis functions. This approximation scheme yields a suboptimal control for the original problem. The optimal control problem is converted to an approximate nonlinear optimization problem with a finite number of decision variables, which can be solved by nonlinear programming techniques. The advantage of piecewise-constant approximation scheme is due to its simplicity and convergence [41].
To solve optimal control problem using the control parameterization method, we approximate the control profile as (see Figure 7) where is a given integer, , are knot points, and , are vectors containing the approximate control values. The knot points satisfy The approximate control can be written as where for a given subinterval , the characteristic function is defined by Note that is a piecewise-constant function with potential discontinuities at the points , . These points are called switching times. Reference [40] surveys the key developments in the control parameterization.
4. Computational Issues
Solving optimal trajectory design problems with direct methods yields a large NLP problem. Some computational issues should be implemented to improve efficiency and robustness. Some of these issues are discussed in this section.
4.1. Scaling
An important issue that arises in the solution of NLP is scaling. It is known that a poorly scaled problem can lead to either extremely slow convergence or divergence [27]. When the relative size of variables or constraints in a problem is vastly different, the variables should be transformed into a relatively similar scale. In trajectory optimization problems the scaling process can be implemented manually by introducing dimensionless distance unit , time unit , velocity unit , and mass unit .
4.2. Sparse Matrix
A sparse matrix is a matrix filled primarily with zeros; the opposite is referred to as dense matrix. One of the major advantages of using sparse matrix is that finite memory can be saved. In the case of direct collocation methods, it is well-known that the vast majority of the Jacobian and Hessian matrix elements are zero. When a large number of nodes are taken to have a smooth state and control time histories, “out-of-memory” problems can be avoided by using sparse matrices. Another advantage is that the computational time is reduced because only nonzero elements need to be manipulated. Some NLP solvers, such as SNOPT [42] and IPOPT [43], take advantage of sparse techniques. Furthermore IPOPT computes second derivatives with quasi-Newton methods.
4.3. Differentiation
Direct methods use gradient-based techniques for solving the NLP, which requires the gradient information of the objective function and constraints with respect to the NLP variables. It is known that providing the analytic derivative results in exact and fast optimization, but sometimes it is impossible or time consuming to compute derivatives analytically. Alternative methods are needed to obtain the necessary gradients. Finite differences are implemented to compute approximate derivatives. Forward and central difference schemes are where is a properly chosen perturbation, which should be sufficiently small to provide a good approximation to the derivatives but not too small to induce round-off errors.
5. Practical Examples
5.1. A Simple Optimal Control Problem
In this section, a simple optimal control problem [44] is solved with Hermite-Simpson method, and certain issues are set to speed-up the computation time and improving accuracy. Consider the problem of finding with fixed initial and final time to minimize the cost subject to system equations and to the boundary conditions with the following bound box of state and control variables The analytical solution of this problem is
This problem is solved by using Matlab’s “fmincon” with “interior-point” algorithm [45]. The analytic solution and final optimized results are plotted in Figures 8(a) and 8(b), where the green line indicates analytic solutions (i.e., and ), and the discrete results obtained from NLP solver are marked by blue asterisk and red cross, respectively. Figure 8(c) shows the error of and in detail, it can be seen that the error of is much lower than that of ( for versus for ). This is because the third-order Hermite interpolation polynomial is implemented to approximate the states, while linear interpolation is used to represent the control curve. In order to check the control accuracy, (66) are integrated from the initial condition with the discrete optimal solution , and then the error of is plotted in Figure 8(d). It can be found that the error of in the whole time domain is below .
| (a) Discrete solution and analytical time history of |
| (b) Discrete solution and analytical time history of |
| (c) Error of and located at discrete points |
(d) Reintegrate dynamics with discrete control solution
The “finite differences” method is the default option to obtain the gradient information, however computational efficiency can be remarkably improved by adding analytical Jacobian and Hessian. The comparison of these three different cases is listed in Table 3. All of the three cases use 100 discrete points. In case 1, the optimization takes 116 steps and is 17.21 s. In case 2, adding Jacobian only does not reduce the number of iterations but saves almost half of the CPU time (17.21 s versus 8.77 s). In case 3, providing analytical Jacobin and Hessian highly reduces both the number of iterations and the computational time, which drops dramatically down to 0.38 s.
The structure of Jacobian is shown in Figure 9(a), where the rows are constraints while the columns are NLP decision variables. The diagonal elements are Jacobian of collocation constraints and the last two lines are Jacobian of initial boundary constraints. Figure 9(b) illustrates the structure of Hessian, where both the row and the column represent NLP decision variables. It can be seen that Hermite-Simpson method has highly sparse Jacobian and Hessian structures.
(a) Jacobian structure, 100 points
(b) Hessian structure, 100 points
Because only few elements in Figures 9(a) and 9(b) are nonzero, sparse matrices can be used. In Figure 8, the accuracy of state and control solution is not so good. A finer mesh with 5000 points is used with “interior-point” method and matrix sparsity, which is impossible with other fmincon algorithms (that use dense matrixes). The results are illustrated in Figure 10. Compared to Figure 8, it is noted that the error of reduces to around . In order to check the accuracy of control solution, the dynamics are reintegrated with the initial condition and optimal control history; the maximal error on drops down to (see Figure 10(d)). A parametric analysis has been carried out with variable grid points. The outcome is reported in Table 4.
| (a) Discrete solution and analytical time history of |
| (b) Discrete solution and analytical time history of |
| (c) Error of and located at discrete points |
(d) Reintegrate dynamics with discrete control solution
5.2. Planar Low-Thrust Orbit Transfer
For a transfer trajectory to be determined, using Cartesian coordinates is the simplest but most disadvantageous choice; this is because a lot of discrete points will be used to catch the rapidly changing position and velocity variables. On the contrary, slowly changing state variables make the NLP problem both efficient and robust. Polar coordinates are then used for the planar two-body dynamics; that is, In (70), is the spacecraft radius, is phase angle, and are the radial and transversal velocities, respectively, is the gravitational constant, is the propulsive acceleration, and is thrust angle (see Figure 11). The initial boundary conditions set the spacecraft in its initial circular orbit at : while the final boundary conditions place the spacecraft into a target circular orbit at : The nondimensional and maximum thrust acceleration are In the following, the high efficient NLP solver IPOPT [43] is used for large-scale nonlinear optimization.
5.2.1. Time-Optimal Problem
The goal of time-optimal problem is to find the functions and that minimize the performance index under the dynamic constraints and the boundary conditions above. This problem is solved with uniformly spaced points with Hermite-Simpson method. (That is, 2401 NLP variables: the states and the controls at the nodes plus the final time.) Figure 12 shows the transfer trajectory; the cyan dashed line is the initial guess (although the convergence basin of direct methods is larger than that of indirect methods, a good initial guess is essential; in time-optimal problem, tangential thrust can be used as a simple but suitable approach to guess an initial solution), the thin blue line denotes the final transfer orbit, and the thick red line shows the thrust arc. Figure 13 shows the time history of state variables. The final transfer time is 55.5. Figure 14 shows the profile of the optimal control variables , ; it can be seen that the engine is on duty along the entire orbit, to shorten the transfer time.
5.2.2. Fuel-Optimal Problem
In the fuel-optimal problem, the functions and are sought to minimize the performance index Figure 15 illustrates the optimal transfer trajectory. This solution has been obtained with mesh points (i.e., 4801 NLP variables). In this problem, tangential thrust with magnitude of is used to produce the initial guess (cyan dashed line in Figure 15). It can be seen that the thruster is on duty across the periapsis; the only maneuver performed at the apoapsis injects the spacecraft into a nearly circular orbit to acquire the final orbit. The final time of flight is 122.3. Figure 16 illustrates the time history of the state variables, whereas Figure 17 illustrates the time history of control variables. A bang-bang structure is found, indicating the optimality of the solution found.
6. Conclusions
In this survey, some issues related to the direct transcription and collocation for optimal low-thrust space trajectory design are discussed. The principles of direct transcription and collocation are given, and emphasis is put on high-order numerical integration. Hermite-Simpson, high-order Gauss-Lobatto, and the pseudospectral schemes are introduced. For the Gauss-Lobatto method, a scheme to handle arbitrary order is given, together with the procedure to derive analytical gradients. The pseudospectral method has an impressive convergence rate; however, it is more suitable for smooth problems. A simple example (two states, one control) has been solved with the Hermite-Simpson method. This example shows outstanding efficiency when analytical gradient is provided. A second example consists in low-thrust transfers in the two-body model. From this example it can be concluded that compared with pseudospectral method, the Gauss-Lobatto method is more suitable for fuel-optimal problems with discontinued control structures. In both examples, effort is put to give insights on the computational issues for a practical implementation of the method. However, if a long duration problem is encountered (with increasing number of on/off structures) or a spacecraft equipped with very small thrust engines is dealt with (large number of spirals), the NLP solver becomes ill-conditioned. In these cases, indirect approach with homotopic or hybrid methods could be more suitable.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
C. Zhang would like to acknowledge the Chinese Scholarship Council (CSC) which financially supported his visit to Politecnico di Milano.