Low Thrust Minimum Time Orbit Transfer Nonlinear Optimization Using Impulse Discretization via the Modified Picard–Chebyshev Method

Low Thrust Minimum Time Orbit Transfer Nonlinear Optimization Using Impulse Discretization via the Modified Picard–Chebyshev Method

Darin Koblick1,2,3,Shujing Xu4,Joshua Fogel5and Praveen Shankar1

The Modified Picard-Chebyshev Method(MPCM)is implemented as an orbit propagation solver for a numerical optimization method that determines minimum time orbit transfer trajectory of a satellite using a series of multiple impulses at intermediate waypoints.The waypoints correspond to instantaneous impulses that are determined using a nonlinear constrained optimization routine,SNOPT with numerical force models for both Two-Body and J2perturbations.It is found that using the MPCM increases run-time performance of the discretized low-thrust optimization method when compared to other sequential numerical solvers,such as Adams-Bashforth-Moulton and Gauss-Jackson 8thorder methods.

LowThrust,Picard–Chebyshev,Optimization,Multi-Impulse,Boundary Value Problem.

1 Introduction

Low thrust propulsion systems provide an efficient means of orbit transfer for both interplanetary trajectories as well as Earth-centric trajectories where mission time is not significantly important[Keaton(2002)].Low thrust propulsion systems,however,operate continuously during the mission and hence depend on having an optimal trajectory defined for the orbit transfer.The optimization of low thrust trajectories using indirect methods encounter highly sensitive two-point boundary value problems with unknown co-state boundary conditions.In the case of direct methods the control history is parameterized leading to a high-dimensional nonlinear programming problem.This results in poor convergence results as the dimensionality of the parameter space grows larger.Thus both direct and indirect methodstypically lead to challenging numerical difficulties.There have been a number of different approaches proposed to solve the general optimization problem of low constant thrust trajectories between orbits[Betts(2000);Kechichian(2007);Yang(2009);Falck and Dankanich(2012);Graham and Rao(2015)].These approaches generally fall into either direct or indirect methods.The direct methods typically convert the low thrust trajectory into a parameter optimization problem which can be solved using a nonlinear programming method.Indirect methods are developed around calculus of variations or the maximum principle[Bryson and Ho(1975)].Using an indirect method can be difficult as there are typically six(resp.four)costate variables for three-dimensional(resp.two-dimensional)cases,which must be determined and they are not only extremely sensitive but difficult to guess for generic orbit transfer problems[Ulybyshev(2009)].An extensive review of the low thrust trajectory optimization based on direct,indirect and stochastic methods and its application can be found in[Kim(2005)].

The continuous low thrust trajectory optimization problem can be discretized into waypoints or nodes with the solution consisting of impulses computed at each of the waypoints using an optimization routine over the entire trajectory.Such a method is described as multi-impulse low thrust trajectory optimization.Similar to the performance of most well behaved numerical methods,the discretized solution approaches the continuous solution as the number of nodes are increased.The method of using multiple impulses to approximate a continuous low thrust trajectory has been addressed by several studies[Xiao-yong,Yi-jun,Hong-bo,and Guo-jian(2013);Polsgrove,Kos,Hopkins,and Crane(2006)].

One of the requirements of the multi-impulse method of low thrust optimization is a numerical solver to propagate the orbit between the waypoints since an analytical solution may not exist for most transfer problems utilizing numerical force models.The choice of this numerical solver plays a significant role in the run-time performance of the multi-impulse method,especially when the number of waypoints are increased to improve accuracy.This paper presents the implementation of the Modified Picard–Chebyshev Method as the numerical solver for orbit propagation within a multi-impulse method of low thrust optimization.The low thrust optimization utilizes the solution from a direct or indirect optimization[Kechichian(1997a)]method as an initial guess and discretizes this guess into evenly distributed waypoints.The method then applies a nonlinear optimization routine,the S-parse Nonlinear OPTimizer(SNOPT)[Gill,Murray,and Saunders(2002)],to determine the instantaneous change in velocity,∆V,at each of these waypoints while minimizing the total∆V over the entire trajectory.It is important to emphasize here that the novelty in the proposed method is the utilization of the Modified Picard–Chebyshev Method[Bai(2010)]for both the boundary value and initial value problem solutions required for the optimization.This implementation im-proves the run-time performance as compared to typical numerical solvers,such as the Adams–Bashforth–Moulton(ODE113)[Bradley(2012)]and Gauss–Jackson 8th(GJ8)[Berry and Healy(2004)]methods.

The rest of the paper is organized as follows:Section 2 addresses some of the fundamental concepts and background.Section 3 describes the multi-impulse low thrust optimization method that utilizes the MPCM as a numerical solver for the initial value and boundary value problems.Section 4 presents the results of utilizing the method on two orbit transfer problems.The results are divided into two parts.The multi-impulse method is compared to an indirect optimization method proposed by[Kechichian(1997a)]and a numerical method based on Edelbaum’s analytical circular transfer problem[Kechichian(1997b)]to show the fidelity of the discretization strategy with the implementation of the Modified–Picard Chebyshev Method.Then,the run-time performance of the method with the Adams-Bashforth-Moulton and Gauss–Jackson 8thorder solvers is compared to that of the implementation with the MPCM.Section 5 provides a discussion of the results and the characteristics of the method and is followed by conclusions and future work.

2 Background

This section discusses some of the fundamental concepts that will be utilized in developing and implementing the multi-impulse low thrust orbit transfer method.

2.1 Optimal Control Problem

Consider the nonlinear system

The optimal control problem is defined as the process of determining the control input u∗such that an associated cost function

is minimized subject to the constraints

In the case of spacecraft trajectory control,these constraints are generally related to the maximum thrust available or the total time of the mission.This optimal control problem can be solved by both direct and indirect methods.In direct methods,the optimal control problem is converted to a nonlinear programming problem which can be solved using a variety of methods.While the direct methods tend to be highly robust,the solutions obtained are generally not very accurate[Von Stryk and Bulirsch(1992)].

Indirect methods are based on calculus of variations and Pontryagin’s minimum principle[Lewis,Vrabie,and Syrmos(2012);Geering(2007)].These methods modify the cost function to include the equations of motion and state and control variable constraints.The Hamiltonian function,H,that results from applying the calculus of variations to the augmented cost function under the constraints can be written as

where λ is a vector of costates associated with the equation of motion constraints.The optimality conditions are obtained by computing the first variation of J resulting in

The above equations represent a Boundary Value Problem(BVP)and while the boundary conditions of the state variables are given, the boundary conditions for the costates are normally unknown.Optimization does not require minimization of the costates λ,but their determination is required for satisfying optimality constraints through Eq.4.There are several methods that have been discussed for solving such BVPs,for example,the shooting method;but since the initial values of the costates are unknown,these methods become extremely cumbersome to implement.

2.2 Orbit Equations of Motion

Equinoctial orbit elements are frequently used in solving low thrust optimization problems as they help in avoiding singularities involved with using the classical Keplerian orbital elements[Chobotov(1996);Broucke and Cefola(1972)].The mapping to convert classical elements to equinoctial elements is

The two-body equations of motion for low thrust based on the equinoctial orbital elements can then be written as

where

The Euler–Lagrange equations can then be determined through numerical integration of

Note that the partial matricies,are provided in the appendicies of[Kechichian(1997a)].

2.3 Low Thrust Optimal Orbit Transfer Problem

The low thrust orbit transfer method can be described as a mission that continuously utilizes low thrust to transfer a satellite from an initial position to a final position.The unique nature of such a mission as compared to an impulse transfer is that the orbit(orbital elements)continuously changes during the mission.The constant low thrust optimal control problem involves the determination of the unit vectorˆu along the spacecraft trajectory that minimizes either transfer time or total∆V with the assumption that the thrustis held constant over the transfer.

This section reviews two methods for constant low thrust transfer optimization based on Euler–Lagrange multipliers and Edelbaum’s method.These two methods are also used as a benchmark for the proposed multi-impulse method.

2.3.1 Minimum Time Continuous Thrust Using Euler–Lagrange Multipliers

Under certain assumptions,there exists a function of adjoint variables λ(t)and a vector function v(t)in which the first order necessary conditions for optimality are valid;see[Von Stryk and Bulirsch(1992)].This is known as the Hamiltonian function of the multi-point boundary value problem defined by Eqs.12–14 and is shown generically as

where f and g are some force functions which describe a dynamic system of nonlinear ODEs(e.g.f(x(t),u(t),t),and g(x(t),u(t),t)).The corresponding coupled differential equations for t0≤t≤tfare defined as

For the low thrust minimum time,Eq.12 can be expressed as[Kechichian(1997a)]

To minimize the total flight time,and indirectly the∆V,the performance index J can be minimized:

For fixedt0,minimization of tfor maximization of−tfandgives rise to the transversality condition,H(tf)=1[Kechichian(1997a);Kechichian(1997c);Chobotov(1996)]resulting in time-optimality.

2.3.2 Edelbaum’s Solution For Non-Coplanar Circular Orbit Transfers

Edelbaum (1961) developed an analytical solution for the optimization of low thrust orbit transfer problems that could be applied to circular trajectories.This section details a numerical solution based on Edelbaum’s method that is valid for inclination changes of less than 114.6°and assumes that the acceleration magnitude is small enough such that the transfer orbits are circular[Kechichian(1997b);Eagle(2013)].Much like the acceleration unit vector that can be determined by the costates in Eq.11,the acceleration can be found for Edelbaum’s solution as well.In this case,the acceleration unit vector is determined from

Note that the acceleration unit vector shown in Eq.16 is in the equinoctial frameand can be transformed to the inertial frame through the following matrices[Chobotov (1996)]:

where

The initial yaw angle of the acceleration vector β0is defined by

where V0is the magnitude of velocity in the initial circular orbit,and Vfis the velocity in the desired circular orbit(e.g., final orbit after the maneuver).The yaw angle and inclination with respect to time can be found by

The total change in velocity required for the transfer,∆V,is given by

The change in inclination with respect to time is

The instantaneous velocity can be computed at anytime during the orbital transfer

Using the orbital-energy-invariance law(vis-viva equation),the semi-major axis with respect to time can be determined from the velocity magnitude

whereµis the standard gravitational parameter of the central body.

The out of plane yaw-angle,β,and the in-plane pitch angle,ψ are related to the equinoctial frame in the following way:

2.4 Overview of Numerical Force Models for Orbit Propagation

One benefit of using numerical techniques when modeling satellite orbit dynamicsis that it is easy to incorporate additional acceleration components into the equations of motion. Cowell’s formulation for the two-body equation of motion [Vallado and McClain (2001)] is

2.4.1 Modeling J2Perturbations

The second zonal harmonic,J2,is among the most dominant of the harmonics and the perturbation acceleration due to J2is significant for orbits at and below GEO altitude[Larson and Wertz(1999)].Eq.31 shows the magnitude of this term when other zonal harmonics are three orders of magnitude smaller(e.g.,10−6or less):

The acceleration acting on the satellite due to J2effects used in this study is provided with respect to the inertial reference frame in Cartesian coordinates[Schaub and Junkins(2003)]

2.5 Lambert Two-Body Problem

Lambert’s problem involves the solution of the orbit between two position vectors and the time of flight required to propagate a satellite from one position to another.The problem can be formulated in the following way:Givenand t2, findandNote thatis the required velocity for the satellite ato intercepatt2.The velocity required to rendezvous withis described asand therequired is shown asandrespectively.

See Fig.1 for an illustration.is the velocity of the satellite at positionbefore the maneuver andis the velocity of the satellite at position→after the maneuver.

Figure 1:Illustration of a Two-Impulse Orbit Transfer Using a Two-Body Lambert Solver

There are many different methods to solve a Lambert Problem:Gauss’Method[Gauss(1857)],Thorne’s Solution [Thorne (1990)], Battin’s Solution [Battin (1999)],Godal’s Method[Godal(1961)],Lancaster’s Method[Lancaster and Blanchard(1969)],and Sun’s Method [Sun (1979)]. The Superior Lambert Solver[Der(2011)]is used in this analysis to obtain the initial thrust vector direction and magnitude.

2.6 Modified Picard–Chebyshev Method

The vectorized matrix formulation of the Picard–Chebyshev Method is a modern concept. It was originally conceived in Feagin and Nacozy (1983); and later it was adapted by Fukushima to run on a vector computer to solve astrodynamical IVPs and BVPs of first order ODEs [Fukushima (1997); Fukushima (1999)]. Since his findings, it has been evaluated as a viable candidate for parallel orbit propagation[Neta (1998)] and further enhanced through the work of [Bai (2010)]. In this section we provide a brief background on the MPCM, which has been vectorized and parallelized for IVPs on multiple processors in our previous studies [Koblick(2012); Koblick, Poole, and Shankar (2012); Koblick and Shankar (2014)].

The formulation for the Modified Picard–Chebyshev Method presented in this section is applicable to certain types of first order BVPs and IVPs,based on that of

Bai(2010),Bai and Junkins(2011a)and Bai and Junkins(2012a).The MPCM can be used to augment the Lambert solution in solving a two-impulse transfer under additional perturbations,as well as to replace the numerical ODE solver which is needed for propagating the solution from one node to the next.

Consider a system of coupled nonlinear first order differential equations which may represent either an IVP or a BVP:

It is desired to obtain a solutionx(t)to Eq.33.Let N be the number of data points within the interval[t0,tf]at which the solution is calculated. The independent variable,time t,can be transformed to another domain,τ,where τ∈ [−1,1]via

The simple relationship between the two domains iswhere ω2=andAfter substitution,Eq.33 can be written as

For the IVP,the initial conditions f(t0,x)are related to g(τ0,x)through

The Picard iteration which serves to update the initial approximation can be computed through

where

Here′′denotes that the first and last terms of the summation are divided by two.For the BVP where the initial and final positions are known,the coefficients needed to find the updated positions with respect to time can be computed from

To compute the updated velocity components(the first derivative of position with respect to time),different coefficients are needed to approximate the Picard iteration:

where

The Picard iteration is performed for j iterations until the relative difference between solutions is less than some desired error tolerance ε,i.e.,

3 Optimal Multi-Impulse Constant Low Thrust Orbit Transfer Method

Like any unperturbed trajectory,a continuous low thrust trajectory can be approximated by discretization.For higher continuous acceleration levels,the discretization can be chosen with respect to the true anomaly while for lower continuous acceleration levels,the discretization can be chosen with respect to time.The multi-impulse discretization approach can be simply described as approximating a continuous thrust trajectory with multiple impulses.Similar techniques have been described and suggested by Xiao-yong,Yi-jun,Hong-bo,and Guo-jian(2013)and Polsgrove,Kos,Hopkins,andCrane(2006).The novelty in this proposed method is to use the Modified Picard–Chebyshev Method not only as a BVP solver,but also an IVP solver within the nonlinear constraint function to support multiple gravitational bodies and complicated numerical force models.The steps necessary in the optimization of a continuous low thrust problem when using the multi-impulse method are described below:

·Step One:Low-Thrust Orbit Transfer ApproximationEdelbaum’s solution for Non-Coplanar Circular Orbit Transfers is used as a good initial approximation for solving continuous low thrust trajectories.It is important to note that there are limitations on the magnitude of the acceleration since the analytic solution starts to break down as the orbit transfer becomes more and more elliptic.In certain cases better initial solutions available through indirect methods such as Kechichian’s Euler–Lagrange Multipliers are used.

·Step Two:Discretization of Orbit Transfer ApproximationThe orbit transfer trajectory developed in step one is discretized into N nodes.Either a two-body boundary value problem solver,such as Lambert’s solution,or another numerical boundary value problem solver such as the Modified Picard–Chebyshev Method is used to compute the instantaneous impulse required to travel from point rN−1to point rN.If the MPCM is used,the Lambert solver is used as an initial guess for the MPCM.

·Step Three:Nonlinear OptimizationGiven the desired constraints developed in the Constraints function,a nonlinear optimization routine,such as SNOPT[Gill,Murray,and Saunders(2002)],fmincon[Waltz,Morales,Nocedal,and Orban(2006)],or IPOPT[Wächter and Biegler(2006)],is used to determine the∆V necessary for each impulse,as well as the time of flight in order to minimize the objective function(which can be either time of flight or total∆V).A numerical solver(MPCM,Adams–Bashforth–Moulton or Gauss-Jackson 8thorder)is utilized to propagate the orbit between waypoints within the optimization process.

·Step Four:Solution Re finementOnce a coarse solution is obtained from step three,this sequence can be repeated using the coarse solution as the initial guess.

Since this approach takes a continuous trajectory and discretizes it into multiple segments,it creates a perfect application for the Modified Picard–Chebyshev Method(the MPCM performs best on partial orbital segments as opposed to many revolutions[Koblick and Shankar(2014)]).The primary result was that the runtime performance of SNOPT could be significantly improved when using the MPCM for numerical propagation as compared to other solvers.

3.1 Initial Set-up

The time of flight,tf,is discretized into N intervals.The impulse vector for each interval∆Viis represented aswhere φ and δ are spherical azimuth and elevation angles with respect to the ECI coordinate frame andηiis a multiplier which should approach one when the acceleration is evenly distributed around the orbit transfer.The optimization variables are

Using a spherical coordinate system allows the search space to be within the range of 0–π radians in elevation and 0–2π radians in azimuth.Upper and lower bounds on ηican be provided to remove very large impulses out of the feasible region.

3.2 Objective Function

This method can be utilized to minimize either the total∆V or the time of flight tf(but not both as the maximal∆V corresponds to the minimal tffor the multiple impulse approximation).The performance index for the minimum∆V objective is

Subsequently,the minimum time of flight performance index is described in Eq.17.For the optimization performed in the following examples,the minimum time of flight was chosen as the objective function with a penalty for∆Viwhen larger than the maximum ∆Vifor the segment,i.e.,when ηi＞ 1.

3.3 Constraints

The constraints associated with an orbital maneuver depend on which objectives are important to the designer.The initial conditions will remain the same regardless of which optimization objectives one has in mind:

For full orbital rendezvous optimization,it may be necessary to fix the final state vector such that the position and velocity at the end of the transfer match the target spacecraft,

For orbit transfers similar to those described by Kechichian,it may be desirable to fix only five of the six orbital elements as show in Eq.51(this can also be written in equinoctial elements to avoid singularities):

For both the minimum∆V and minimum tfoptimization problems,the low thrust constraint is satisfied by

where ηiis the maximum ratio(typically around 1)where

3.4 Optimization Routine

A 64-bit MATLAB executable of SNOPT,was used for the minimum time of flight optimization shown in the results section.While both SNOPT and fmincon were used to solve the same constrained nonlinear optimization problems,SNOPT showed better convergence characteristics for this specific implementation.

4 Results

The optimal method developed in this paper is used to perform low thrust orbit maneuvers for two cases:

Case 1Minimum Time Non-Coplanar Circular Transfer

Case 2Minimum Time LEO to GEO Transfer

4.1 Non-Coplanar Circular Transfer

For this example, the constant low thrust transfer of a satellite from inclined circular GEO to equatorial circular GEO is optimized.Initial conditions are a0=42000 km,with a constant accelerationThe final conditions of the target orbit are af=42000 km,ef=0.00,if=0.00°;Ωf,ωf,and Mfare determined by the initial approximation;tfis to be minimized.Using Eq.24,the analytical solution to this problem corresponds to a∆V of 0.8419 km/s and a time of flight of 84193.88 seconds.The effects on optimization when a different number of impulses are chosen with respect to each orbit revolution are shown in Table 1.It is interesting to note that there is an optimal setting to achieve high feasibility and the most realistic low thrust segmentation(η≤1.000);this occurs around eight impulses/revolutions.It should be noted that this paper uses the same de finitions for optimality and feasibility as outlined in the SNOPT User’s Guide[Gill,Murray,and Saunders(2006)].

Table 1:Number of Impulses vs Total∆V,Time of Flight,and Feasibility.

A top view of the noncoplanar orbit transfer is shown in Fig.2;Spline interpolated thrust vectors are overlaid in Fig.3.

Figure 2:Birds Eye View of Transfer

Figure 3:Acceleration Vectors

4.1.1 Comparison to Numerical Solution Based on Edelbaum’s Method

The solution obtained by the multi-impulse optimization method is compared to an implementation of Edelbaum’s method and a trial copy of GPOPS-II MATLAB software.GPOPS-II is general-purpose software designed to solve multiple phase optimal control problems using variable-order Gaussian quadrature collocation methods [Patterson and Rao (2014)]. The low-thrust orbit transfer example that was provided with a copy of GPOPS-II(developed from a published optimization example[Betts(2001)])was Modified to minimize time of flight instead of fuel.Table 2 shows the results for the Edelbaum Transfer,the Multi-Impulse Transfer,and the GPOPS-II transfer case.In this case,the optimization routine minimizesthe total time of flight required for the transfer instead of the∆V.Theoretically,the minimum transfer time should correspond to the minimum∆V for the low continuous thrust case.Multiplying the total transfer time by the acceleration should produce the∆V associated with the transfer.As observed from the table,the solution will approach that of the constant acceleration case as the total number of impulses increase.

Table 2:Final Orbital Elements and Flight Times for GEO-GEO Transfer.

4.1.2 Numerical Solution Using Additional Perturbations

For the J2perturbation force model as described in Eq.32,the Multi-Impulse Transfer solution is obtained by running SNOPT,shown in Table 3.Table 4 shows the run-time performance results associated with different numerical propagation techniques within the nonlinear constraint routine.It can be seen that the implementation with MPCM results in a significant improvement in performance for both the two-body and J2perturbations as compared to the implementations with Adams-Bashforth-Moulton and Gauss-Jackson 8thorder methods.In Figs.4–6,the thrust vectors associated with the transfer for the two-body solution are shown.Note that with a minor perturbation such as J2,these vectors may show large changes in direction,suggesting that this optimization problem is not well posed.

4.2 LEO to GEO Minimum Time Low Thrust Transfer

This section addresses the solution of a LEO to GEO transfer scenario based on Kechichian[1997a].The initial conditions for this scenario are a0=7000 km,e0=0.00,i0=28.5°,Ω0=0.00°,ω0=0.00°,M0= −220°,with a constant accelerationAnd the final conditions of the target GEO are af=42000 km,ef=0.001,if=1.00°,Ωf=0.00°,ωf=0.00°,Mf= −43.779715°,and tf=58624.094 sec.This equates toNote that Kechichian allowed the final mean anomaly to be free in order to optimize where in the GEO orbit the satellite may enter.In this paper,it was chosen to fix this entry point at his optimal two-body solution of 43.7797°.This allows for a full rendezvous from one fixed inertial position to another while minimizing time of flight and subsequently∆V,assuming that the final implementation of thrust will be constant and continuous.A comparison of convergence performance versus total number of impulses is provided in Table 5.While fewer impulses may appearto be attractive for this optimization procedure,the resulting trajectory is highly eccentric and is not easily converted to a continuous thrust solution.30 impulses were chosen for the optimization process as it meets the required feasibility and accuracy.

Table 3:Final Orbital Elements and Flight Times for GEO-GEO Transfer with J2 Effects.

Table 4:SNOPT Run Time Performance of GJ8,ODE113,and MPCM(GEOGEO).

Table 5:Number of Impulses vs Total∆V,Time of Flight,and Feasibility.

Figure 4:X-axis

Figure 5:Y-axis

Figure 6:Z-axis

4.2.1 Comparison to Kechichian Solution

With the Lagrange multipliers specified by Kechichian as input into the low thrust orbit propagation IVP(see Table 6),this trajectory was discretized into 30 nodes which were then optimized through SNOPT.The∆V was computed as 5.743 km/s.The minimized flight time was reported as 58624.093 seconds.Table 7 compares the final position of the satellite after the low thrust transfer.Fig.7 depicts a three dimensional view of the full orbit transfer.

Table 6:Initial Costate Values Corresponding to the LEO-GEO Transfer Problem Outlined by Kechichian[1997a]and reprinted in Chobotov[1996].

Table 7:Comparison of Final Orbital Elements.

Figure 7:30-Node ECI Thrust Vector Comparison Multiple Impulse vs Costates(Blue:Multi-Impulse Method;Black:Costate Method)

The unit vectors of the acceleration in the ECI reference frame are compared at each node along the transfer trajectory in Figs.8–10.Note that the black dot represents the solution found from optimizing the costates while the blue line is the thrust vector at each discrete point as found from running the multiple impulse method.

Figure 8:x-axis

Figure 9:y-axis

4.2.2 Numerical Solution Using Additional Perturbations

For the J2perturbation force model as described in Eq.32,Table 8 shows the Multi-Impulse Transfer solution obtained by running SNOPT.

Table 8:Final Orbital Elements and Flight Times for LEO-GEO Transfer with J2 Effects.

Table 9:SNOPT Run Time Performance of GJ8,ODE113,and MPCM(LEOGEO).

Table 9 shows the comparison of the run times of the multi-impulse low thrust optimizations with implementation of MPCM,ODE113 and GJ8.The results show a significant improvement in the performance of the MPCM as compared to the other two solvers.

5 Discussion

This study presents the computational run-time performance gain of using the MPCM to solve discretized constant low thrust orbit transfer optimization problems in conjunction with a sparse sequential quadratic programming(SQP).Previous work in relation to the MPCM has shown applications for Lambert targeting problems,two-body coplanar low thrust trajectory optimization,and Station-Keeping of Translunar Halo Orbits[Bai and Junkins(2011b);Bai and Junkins(2012b)].

In regard to performance on numerical orbit propagation problems,previous research has shown that when distributed over many GPUs,the MPCM can reduce the computational time between one and two orders of magnitude over other conventional sequential ODE solvers such as Runge–Kutta 4(5)and Runge–Kutta–Nystrom 12(10)[Bai and Junkins(2012a)].Additionally,the work by Koblick and Shankar(2014)has shown that a parallel MPCM implementation across just five CPU cores can reduce the computational time of orbit propagation to below the performance of the fastest sequential numerical ODE solvers(e.g.,Runge Kutta 4(5),Gauss Jackson 8,and Dormand and Prince(8,7))[Koblick and Shankar(2014)].Fukushima,after testing a vectorized version of the Picard–Chebyshev method on a vectorized computer,theorized that the performance improvement of solving perturbed orbit propagation problems can be between two and three orders of magnitude[Fukushima(1999)].This suggests it is likely that the run-time performance of this low thrust discretization method,when using the MPCM,can benefit significantly when switching to a parallel implementation(especially when the force model is increased in its complexity).

Direct optimization,for moderate dimensional parametric representation of the unknown control function,is widely known to efficiently lead to good sub-optimal results,and the dimensionality can be frequently increased to re fine the solution without becoming numerically intractable.Excellent approximations of the optimal trajectory can be obtained by using a direct approach,and this approximate truth is con firmed by the results obtained in this paper.

If we upgrade our numerical force models to include multi-body perturbations,this will allow the multi-impulse method to be extended to interplanetary missions where,traditionally,patched conics and multi-segmented optimization would need to be performed.There are some unique characteristics to consider when using a direct optimization method which discretizes the problem in order to optimize an orbit transfer:

·Very good initial guess of trajectory required–because there are many segments,it may take a large number of function evaluations(if feasible)for a nonlinear solver to converge to an optimal transfer as the search space is significantly large(N×3+1)dimensions.

·Segments must be a fraction of each orbit revolution–as shown in the thrust restriction identified in Eq.52,the fewer number of segments used,the faster the optimization(less optimization variables reduce the dimensionality of the problem),but the solution will approach a high thrust transfer and will no longer be similar to that of a low thrust trajectory.

·Real time trajectory controller extension–with a finite number of impulses,the solution is not directly translatable to instantaneous thrust vectors with respect to time.The trajectory must be modeled with some type of trajectory controller,such as a predictive method.

6 Conclusions and Future Work

This paper presents a multi-impulse low thrust trajectory optimization method that utilizes the Modified Picard–Chebyshev Method to improve the run-time performance of the nonlinear optimization routine even with the inclusion of perturbation terms in the gravitation force model.The method shows a significant improvement in the run-time performance as compared to ODE113 and GJ8.The multi-impulse method with the implementation of MPCM was also compared to well known solutions and tools indicating similar accuracies.Future efforts will be directed towards parallelizing the Modified Picard–Chebyshev Method when called from a nonlinear optimization solver in order to further increase the performance benefits of using a perturbative approach to solving numerical propagation problems.This will be followed by an evaluation of existing approaches to convert a multi-segmented solution into a trajectory control algorithm which may be used on-board a spacecraft once an optimal transfer has been found.This multi-impulse discretization method could also be applied to the N-body problem,specifically for interplanetary low thrust rendezvous,low thrust gravity assists(LTGA),or multiple-LTGA problems.One approach used to solve the rendezvous and LTGA problems exploits the natural dynamics of N-body models in the formation of an initial guess trajectory.LTGA consist of a set of patched controlled planar Circularly Restricted N-Body Problems(CRNBP)which analytically generate low energy initial guess formulations(e.g.,invariant manifolds).This initial guess is then fed to a traditional N-body optimizer to compute a complete trajectory,inclusive of planetary inclinations and eccentricities.This approach could extend the multi-impulse discretization method by using the CRNBP as an initial approximation in place of the Edelbaum or the Euler–Lagrange indirect method.Combining systems of CRNBPs as the initial guess to the multi-impulse discretization method could yield a robust strategy for computing complex interplanetary transfer trajectories.

Acknowledgement:The authors would like to thank Phillip Gill from the U-niversity of California,San Diego for providing a copy of the Sparse Nonlinear OPTimizer(SNOPT)software which was used during this study.They would also like to thank Anil Rao from the University of Florida for providing a trial copy of GPOPS-II which was used for comparison results.Also,the authors thank Professor John Junkins from Texas A&M University for his review and feedback of this paper.

Bai,X.(2010): Modified Chebyshev-Picard Iteration Methods For Solution Of Initial Value And Boundary Value Problems.PhD thesis,Texas A&M University,2010.

Bai,X.;Junkins,J.(2011a): Modified chebyshev-picard iteration methods for orbit propagation.The Journal of the Astronautical Sciences,vol.58,no.4,pp.583–613.

Bai,X.;Junkins,J.(2011b): Modified chebyshev-picard iteration methods for solution of boundary value problems.Journal of the Astronautical Sciences,vol.58,no.4,pp.615–642.

Bai,X.;Junkins,J.(2012a): Modified chebyshev-picard iteration methods for solution of initial value problems.Journal of the Astronautical Sciences,vol.59,no.1,pp.327–351.

Bai,X.;Junkins,J.(2012b): Modified chebyshev-picard iteration methods for station-keeping of translunar halo orbits.Mathematical Problems in Engineering,vol.2012,pp.1–18.

Battin,R.(1999): An Introduction To The Mathematics And Methods Of Astro-dynamics.AIAA Education Series.

Berry,M.;Healy,L.(2004): Implementation of Gauss-Jackson integration for orbit propagation.Journal of the Astronautical Sciences,vol.52,pp.331–357.

Betts,J.(2000):Very low-thrust trajectory optimization using a direct sqp method.Journal of Computational and Applied Mathematics,vol.120,no.1,pp.27–40.

Betts,J.(2001): Practical Methods For Optimal Control Using Nonlinear Programming.SIAM, first edition.

Bradley,B.(2012): Numerical integration in astrodynamics(initial value problem). http://ccar.colorado.edu/asen5519/imd/documents/Integration_Lecture_Slides.pdf,2012.

Broucke,R.;Cefola,P.(1972): On the equinoctial orbit elements.Celestial Mechanics,vol.5,no.3,pp.303–310.

Bryson,A.;Ho,Y.-C.(1975): Applied Optimal Control:Optimization,Estimation and Control.CRC Press, first edition.

Chobotov,V.(1996): Orbital Mechanics.AIAA Education Series,second edition.

Der,G.J.(2011): The superior lambert algorithm.In Advanced Maui Optical and Space Surveillance Technologies Conference,Maui,HI.

Eagle,D.(2013):Low thrust transfer between non-coplanar circular orbits.Online,2013.

Edelbaum,T.(1961):Propulsion requirements for controllable satellites.American Rocket Society Journal,vol.31,pp.1079–1089.

Falck,R.;Dankanich,J.(2012):Optimization of low-thrust spiral trajectories by collocation.In Astrodynamics Specialist Conference,no.2012-4423.AIAA/AAS.

Feagin,T.;Nacozy,P.(1983): Matrix formulation of the picard method for parallel computation.Celestial mechanics,vol.29,no.2,pp.107–115.

Fukushima,T.(1997): Vector integration of dynamical motions by the picard-chebyshev method.The Astronomical Journal,vol.113,no.6,pp.2325–2328.

Fukushima,T.(1999):Parallel/vector integration methods for dynamical astronomy.Celestial Mechanics and Dynamical Astronomy,vol.73,no.1,pp.231–241.

Gauss,K.(1857): Theoria Motus(in Latin)Theory of Motion of the Heavenly Bodies about Sun in a Conic Section.Little,Brown and Company.

Geering,H.(2007): Optimal Control With Engineering Applications.Springer,first edition.

Gill,P.;Murray,W.;Saunders,M.(2002):Snopt:An sqp algorithm for largescale constrained optimization.Journal on Optimization,vol.12,no.4,pp.979–1006.

Gill,P.;Murray,W.;Saunders,M.(2006): Users guide for snopt version 7:Software for large-scale nonlinear programming,2006.

Godal,T.(1961):Method for determining the initial velocity vector corresponding to a given time of free flight transfer between given points in a simple gravitational field.Astronautik,vol.2,pp.183–186.

Graham,K.;Rao,A.(2015):Minimum-time trajectory optimization of multiple revolution low-thrust earth-orbit transfers.Journal of Spacecraft and Rockets,vol.52,no.3,pp.711–727.

Keaton,P.(2002): Low-thrust rocket trajectories.Technical Report LA-10625-MS,Los Alamos National Laboratory,Los Alamos,NM,2002.

Kechichian,J.(1997a):Optimal low-earth-orbit-geostationary-earth-orbit intermediate acceleration orbit transfer.Journal of Guidance,Control,and Dynamics,vol.20,no.4,pp.803–811.

Kechichian,J.(1997b):Reformulation of edelbaum’s low-thrust transfer problem using optimal control theory.Journal of Guidance,Control,and Dynamics,vol.20,no.5,pp.988–994.

Kechichian,J.(1997c):Mechanics of trajectory optimization using nonsingular variational equations in polar coordinates. Journal of Guidance,Control,and Dynamics,vol.20,no.4,pp.812–818.

Kechichian,J.(2007): The streamlined and complete set of the nonsingular j2-perturbed dynamic and adjoint equations for trajectory optimization in terms of eccentric longitude.Journal of the Astronautical Sciences,vol.55,no.3,pp.325–348.

Kim,M.(2005): Continuous Low-Thrust Trajectory Optimization:Techniques and Application.PhD thesis,Virginia Polytechnic Institute and State University,2005.

Koblick,D.(2012): Parallel high-precision orbit propagation using the modified picard-chebyshev method.Master’s thesis,California State University,Long Beach,2012.

Koblick,D.;Poole,M.;Shankar,P.(2012):Parallel high-precision orbit propagation using the Modified picard-chebyshev method.In International Mechanical Engineering Congress and Exposition,pp.587–605,Houston,TX.ASME.

Koblick,D.;Shankar,P.(2014): Evaluation of the Modified picard-chebyshev method for high-precision orbit propagation.Journal of Aerospace Engineering,vol.28,no.5.

Lancaster,E.;Blanchard,R.(1969): A unified form of lambert’s theorem.Technical Report TN D-5368,NASA,1969.

Larson,W.;Wertz,J.(1999): Space Mission Analysis And Design.Space Technology Library.Microcosm,third edition.

Lewis,F.;Vrabie,D.;Syrmos,V.(2012): Optimal Control.Wiley,Hoboken,NJ,third edition.

Neta,B.(1998):Parallel version of special perturbations orbit propagator.Technical report,Naval Postgraduate School,Monterey,CA,1998.

Patterson,M.;Rao,A.(2014):Gpops-ii:A matlab software for solving multiplephase optimal control problems using hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming.Transactions on Mathematical Software(TOMS),vol.41,no.1,pp.1–39.

Polsgrove,T.;Kos,L.;Hopkins,R.;Crane,T.(2006):Comparison of performance predictions for new low-thrust trajectory tools.In Astrodynamics Specialist Conference and Exhibit,Keystone,CO.AIAA/AAS.

Schaub,H.;Junkins,J.(2003): Analytical Mechanics of Space Systems.AIAA Education Series,second edition.

Sun,F.-T.(1979):On the minimum time trajectory and multiple solutions of lambert’s problem.In American Institute of Aeronautics and Astronautics Conference,Provincetown,MA.AIAA.

Thorne,J.(1990): Series reversion/inversion of lambert’s time function.In Astrodynamics Conference.AIAA.

Ulybyshev,Y.(2009):Spacecraft trajectory optimization based on discrete sets of pseudoimpulses.Journal of Guidance,Control,and Dynamics,vol.32,no.4,pp.1209–1217.

Vallado,D.;McClain,W.(2001): Fundamentals of Astrodynamics and Applications.Space Technology Library.Springer,second edition.

Von Stryk,O.;Bulirsch,R.(1992): Direct and indirect methods for trajectory optimization.Annals of operations research,vol.37,no.1,pp.357–373.

Wächter,A.;Biegler,L.(2006): On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.Mathematical programming,vol.106,no.1,pp.25–57.

Waltz,R.;Morales,J.-L.;Nocedal,J.;Orban,D.(2006): An interior algorithm for nonlinear optimization that combines line search and trust region steps.Mathematical Programming,vol.107,no.3,pp.391–408.

Xiao-yong,J.;Yi-jun,L.;Hong-bo,Z.;Guo-jian,T.(2013):A multi-impulse extended method for low-thrust trajectory optimization.In Recent Advances in Space Technologies(RAST),2013 6th International Conference on,pp.365–368.IEEE.

Yang,G.(2009): Direct optimization of low-thrust many-revolution earth-orbit transfers.Chinese Journal of Aeronautics,vol.22,no.4,pp.426–433.

1California State University-Long Beach,CA,USA.

2Millennium Space Systems Inc.,El Segundo,CA,USA.

3Claremont Graduate University,CA,USA.

4University of California-San Diego,CA,USA.

5University of Southern California,Los Angeles,CA,USA.