Indirect Optimization of Satellite Deployment into a Highly Elliptic Orbit

Abstract

The analysis of the optimal strategies for the deployment of a spacecraft into a highly elliptic orbit is carried out by means of an indirect optimization procedure, which is based on the theory of optimal control. The orbit peculiarities require that several perturbations are taken into account: an 8 × 8 model of the Earth potential is adopted and gravitational perturbations from Moon and Sun together with solar radiation pressure are considered. A procedure to guarantee convergence and define the optimal switching structure is outlined. Results concerning missions with up to 4.5 revolutions around the Earth are given, and significant features of this kind of deployment are highlighted.

1. Introduction

Space trajectory optimization has gained importance thanks to the development of digital computers and the demanding features of space missions that are currently envisaged. Payload maximization is fundamental to ensure mission feasibility and to reduce costs; sometimes, flight-time minimization is sought to comply with operational requirements. Numerical methods for trajectory optimization can be in general classified into three main groups: indirect methods, direct methods, and evolutionary algorithms. The last group is receiving a great attention, as these algorithms are intrinsically apt to multidisciplinary and multiobjective optimization and in principle are capable of achieving the global optimum in a very large search space. When low-thrust trajectories are considered, evolutionary algorithms typically rely on approximations to maintain the computational time at reasonable values, as a large number of evaluations are required to achieve the solution; for these reasons, they usually provide only an estimate of the optimal solution and a refinement is required to obtain the optimum; direct methods are often used for this purpose. As a consequence, low-thrust trajectories are often dealt with by either direct or indirect methods that typically perform single-objective optimization, attaining a local optimum close to a tentative solution; additional effort is necessary to assure the achievement of the global optimum.

An accurate comparison of direct and indirect methods is found in Betts [1]. A wide number of similarities between the methods is highlighted. Direct methods introduce a parametric representation of the control and/or state variables; the large number of variables, which are required to accurately describe the problem, usually leads to long computational times that can be reduced by taking the matrix sparsity into account. Indirect methods are fast due to the reduced number of variables and may offer higher accuracy and interesting theoretical insight into the problem characteristics. However, three main drawbacks of indirect techniques need to be underlined [1]: analytic expressions for the optimum necessary conditions must be derived, the region of convergence for a root-finding algorithm may be small, and, for problems with path inequalities, it is necessary to guess the sequence of constrained and unconstrained subarcs. It is, however, important to note that also direct methods rely on a tentative solution and may not converge to the optimal solution. However convergence difficulties prevent indirect methods from finding a solution, whereas direct methods find at least a suboptimal solution.

The authors have been using an indirect method for many years and have developed a procedure that mitigates the drawbacks of this approach; in particular the formulation of the optimization problem is made quite simple and attention can be paid to strategies to achieve convergence. In the past two decades the procedure has been tested on different and often difficult problems of spaceflight mechanics and very accurate results have been obtained [2–8]. Recently, the authors had to employ their procedure in a complex application concerning the finite-thrust deployment of a satellite into a highly elliptic orbit (HEO). Dynamical model and problem data were assigned; numerical results were requested for the sake of comparison. HEOs are gaining interest in the scientific community because they represent a less expensive alternative to halo orbits around the Earth-Moon Lagrangian points; in fact, they combine great semimajor axis with large eccentricity, thus providing long stays far from the Earth disturbances at the apogee, where spacecraft velocity is very low. For these reasons they are envisaged as proper operational orbits of spacecraft devoted to deep-space observation. An example is the now-canceled Simbol-X mission [9], which was designed to employ a two-spacecraft formation flying on a HEO to create a new-generation X-ray telescope.

HEOs present low perigee and high apogee; thus, perturbations due to the Earth asphericity and gravitational attraction of Sun and Moon are important and cannot be neglected. In the literature great attention is paid to the behavior of a satellite in such perturbed environment in order to find fuel-saving station-keeping strategies and improve lifetime. Less attention has been paid to the optimization of the deployment maneuver. Typically, these complex problems are faced using direct methods [10, 11]; indirect optimization methods have been widely applied to interplanetary trajectories (due to the simple dynamics), whereas their application to geocentric problems (involving perturbations) is less frequent. In these cases the attention is often focused on continuous-thrust solutions [12–14], and minimum-time transfers are generally sought; minimum-fuel missions with coast arcs would be too long when a high-𝐼_sp very-low-thrust engine is employed. A very-low-thrust multirevolution transfer was analyzed [15] using averaging techniques to deal with perturbations and a continuation scheme to achieve the optimal solution.

Spacecrafts often use chemical engines for orbit deployment. In these cases the thrust level is not low enough for a convenient continuous-thrust transfer and coast arcs are necessary to maintain the fuel consumption within acceptable levels. Several authors [16–18] have dealt with this kind of bang-bang problems but perturbations are usually neglected. In a recent article [19], Thevene and Epenoy include J2 perturbations, whose effects are exploited to reduce the propellant consumption in reconfiguring a four-spacecraft formation; Chuang et al. [20] discuss the effects of atmospheric drag and Earth oblateness for a fixed-duration transfer.

In the present paper the multirevolution finite-thrust deployment of a satellite into an assigned HEO is optimized, taking the relevant perturbations into account. The dynamical model considers an 8 × 8 model of the Earth gravitational potential, solar radiation pressure, and gravitational perturbations from Moon and Sun, whose positions are obtained via ephemerides. The maneuver is essentially a perigee raising and atmospheric drag, which is negligible, can be omitted; however the procedure described here could include aerodynamic forces [4, 21]. The perigee radius exhibits significant variations, even during a ballistic orbit, due to the influence of Sun and (mostly) Moon; these variations are sensitive to the position of the perturbing bodies, thus depending on the departure date. During each apogee passage, the engine and the luni-solar gravitation modify the perigee height; the orbital period is also changed, and this has influence on the position of the perturbing bodies at the following apogees. The effect on the burn structure (i.e., the order of the perigee and apogee burns), which minimizes the fuel consumption, is hard to predict.

The analytic formulation of the indirect optimization problem is quite simple, and the necessary conditions for optimality are derived with little effort. Numerical convergence to the optimal solution is instead difficult to achieve and the dependence of the burn structure on the departure date makes the use of continuation techniques on the departure date unfeasible. The present paper introduces a proper stepped procedure to find out the optimal burn structure for any assigned departure date; numerical examples are presented. A fast convergence to the optimal solution is obtained for a wide range of departure dates; the switching structure of the optimal trajectory highlights how perturbations influence the deployment and how their effects can be exploited to save propellant. Results prove that an indirect method can be effective in analyzing this kind of maneuvers; an accurate problem formulation and a suitable solution procedure help to mitigate the drawbacks of the indirect approach.

2. Dynamic Problem

The paper considers the transfer of a satellite from an elliptic parking orbit with low perigee to the final operational HEO. Reference is made to the now-canceled Simbol-X mission; the characteristics of the initial and final osculating orbits are given in Table 1. The mission starts at the perigee of the initial orbit, where the launcher has released the spacecraft, and ends when the apogee of the final orbit is reached. Initial position, velocity, and mass are known. Only semimajor axis 𝑎, eccentricity 𝑒, and true anomaly 𝜈 are assigned at the final point; no constraints are imposed on the value of the other orbital parameters (deployment to an assigned final orbit was dealt with in [22]). The initial mass is 960 kg, and the final mass is maximized. For operational reasons, thrust cannot be used during the first revolution (orbit acquisition) and in the proximity of the last apogee passage (start of operations).

Table 1 

Initial and final orbit characteristics.

The spacecraft is modeled as a point with variable mass. Position 𝐫, velocity 𝐯,and mass 𝑚 of the spacecraft are the problem state variables, described by differential equations d𝐫d𝑡=𝐯,d𝐯d𝑡=−𝜇𝐫𝑟³+𝐓𝑚+𝐚_𝑝,d𝑚𝑇d𝑡=−𝑐. (1) The trajectory is controlled by the thrust vector 𝐓(the effective exhaust velocity 𝑐 is assumed constant). The perturbing acceleration 𝐚_𝑝=𝐚_𝐽+𝐚_lsg+𝐚_srp (2) is composed of the perturbations due to the Earth asphericity 𝐚_𝐽, luni-solar gravity 𝐚_lsg, and solar radiation pressure 𝐚_srp.

The Earth Mean Equator and Equinox of Epoch J2000 reference frame (i.e., EME2000) is adopted; 𝐈, 𝐉, and 𝐊 are unit vectors along the axes of EME2000. Precession and nutation are neglected. Position is described by radius 𝑟, right ascension 𝜗, declination 𝜑 as 𝐫=𝑟cos𝜗cos𝜑𝐈+𝑟sin𝜗cos𝜑𝐉+𝑟sin𝜑𝐊. (3)

The topocentric reference frame, identified by unit vectors 𝜾(radial), 𝚥 (eastward), and 𝐤 (northward), is introduced. One has 𝜾𝚥𝐤=𝐈𝐉𝐊cos𝜗cos𝜑sin𝜗cos𝜑sin𝜑−sin𝜗cos𝜗0−cos𝜗sin𝜑−sin𝜗sin𝜑cos𝜑. (4) The position vector in the topocentric frame is 𝐫=𝑟𝜾,and the velocity vector is expressed as 𝐯=𝐫̇=𝑢𝜾+𝑣𝚥+𝑤𝐤 (5) with 𝑢, 𝑣, and 𝑤 being radial, eastward, and northward components, respectively. The scalar state equations are easily derived: d𝑟d𝑡=𝑢,d𝜗=𝑣d𝑡,𝑟cos𝜑d𝜑=𝑤d𝑡𝑟,d𝑢𝜇d𝑡=−𝑟²+𝑣²+𝑤²𝑟+𝑇_𝑢𝑚+𝑎_𝑝𝑢,d𝑣=d𝑡−𝑢𝑣+𝑣𝑤tan𝜑𝑟+𝑇_𝑣𝑚+𝑎_𝑝𝑣,d𝑤=d𝑡−𝑢𝑤−𝑣²tan𝜑𝑟+𝑇_𝑤𝑚+𝑎_𝑝𝑤,d𝑚𝑇d𝑡=−𝑐, (6) where subscripts 𝑢, 𝑣, and 𝑤 denote the components along 𝜾, 𝚥, and 𝐤, respectively. It is important to note that the state equations are relatively simple with this set of variables in comparison, for instance, to the use of equinoctial elements. This fact facilitates the analytical derivation of the necessary condition for optimality.

2.1. Earth Potential Model

The Earth potential description is based on the Earth Gravitational Model EGM2008, which provides normalized spherical harmonic coefficients for the Earth gravitational potential; the “Tide Free’’ system is used [23]. The developed code can be quickly modified to consider higher-degree terms or the “Zero Tide’’ system. The Earth’s rotation is assumed to be uniform, neglecting precession and nutation. The EME2000 reference frame is adopted. The gravity model is described in detail in [24].

According to EGM2008, the potential corresponding to the Earth asphericity is expressed as 𝜇Φ=−𝑟^𝑁_𝑛=2𝑟_𝐸𝑟^𝑛𝑛_𝑚=0𝐶_𝑛𝑚cos𝑚𝜆+𝑆_𝑛𝑚sin𝑚𝜆𝑃_𝑛𝑚sin𝜑, (7) where 𝜇 is the Earth gravitational parameter and 𝑟_𝐸 is the semimajor axis of the Earth ellipsoid. In this paper 𝑁 is chosen equal to 8. The associated Legendre functions 𝑃_𝑛𝑚(sin𝜑) and the spherical harmonic coefficients 𝐶_𝑛𝑚 and 𝑆_𝑛𝑚 are used in the unnormalized form that permits faster computations. Normalized quantities would allow for a greater accuracy, which is not necessary for the present application.

The terrestrial latitude coincides with declination 𝜑, as nutation is neglected. The terrestrial longitude 𝜆 is obtained as 𝜆=𝜗−𝜗_𝐺ref−𝜔_𝐸(𝑡−𝑡_ref), where 𝜗_𝐺ref is the Greenwich right ascension at the reference time 𝑡_ref (51544.5 MJD) and 𝜔_𝐸 is evaluated on the basis of the sidereal day, neglecting precession.

The perturbing acceleration due to the Earth asphericity is the gradient of −Φ, and its components in the topocentric frame are thus evaluated as 𝑎_𝐽𝑢=−𝜕Φ,𝑎𝜕𝑟_𝐽𝑣()=−𝜕Φ/𝜕𝜗,𝑎𝑟cos𝜑_𝐽𝑤=−𝜕Φ/𝜕𝜑𝑟. (8) Differentiation with respect to 𝑟 and 𝜗 is straightforward; derivatives with respect to 𝜑 require the derivatives of the associated Legendre functions, which are obtained recursively, exploiting the properties of the Legendre polynomials. Derivatives are evaluated directly with respect to 𝜑 (some authors use the colatitude 𝜋/2−𝜑, the only difference being a sign change of the derivatives); one has, posing 𝑃_𝑛𝑚=0 for 𝑚>𝑛, d𝑃_𝑛𝑚=𝑃d𝜑_𝑛1𝑃for𝑚=0,_𝑛(𝑚+1)−(𝑛+𝑚)(𝑛−𝑚+1)𝑃_{𝑛(𝑚−1)}2for𝑚>0. (9) Further details can be found in [25–27].

2.2. LuniSolar Perturbation

Moon and Sun positions are evaluated using DE405 JPL ephemeris [28], which directly provide the body position in rectangular coordinates 𝑥_𝑏, 𝑦_𝑏, and 𝑧_𝑏 (with either subscript 𝑏=𝑠 for Sun or 𝑏=𝑙 for Moon) with respect to the Earth in the International Celestial Reference Frame and therefore in the EME2000 frame (differences between these frames are very small and can be neglected in the present problem). The perturbing acceleration on the spacecraft, which is caused by a body with gravitational parameter 𝜇_𝑏 and position vector with respect to the Earth 𝐫_𝑏=𝑥_𝑏𝐈+𝑦_𝑏𝐉+𝑧_𝑏𝐊, is given by the difference of the gravitational accelerations that the perturbing body causes on spacecraft and Earth, that is, 𝐚_𝑏𝑔𝜇=−^𝑏𝑅³𝜇𝐑−^𝑏𝑟³_𝑏𝐫_𝑏, (10) where 𝐑=𝐫−𝐫_𝑏 is the spacecraft relative position vector with respect to the perturbing body (and −𝐫_𝑏 is the Earth relative position), as shown in Figure 1.

Figure 1 

Schematic geometry of gravitational perturbations.

The acceleration is projected onto the topocentric frame (based on the spacecraft position) to easily obtain 𝑎_𝑏𝑔𝑢=𝜇_𝑏𝑅³𝑟_𝑏𝑢−𝜇−𝑟^𝑏𝑟³_𝑏𝑟_𝑏𝑢,𝑎_𝑏𝑔𝑣=𝜇_𝑏𝑅³𝑟_𝑏𝑣−𝜇_𝑏𝑟³_𝑏𝑟_𝑏𝑣,𝑎_𝑏𝑔𝑤=𝜇_𝑏𝑅³𝑟_𝑏𝑤−𝜇_𝑏𝑟³_𝑏𝑟_𝑏𝑤 (11) with 𝑅=[𝑟−(𝑟_𝑏)_𝑢]²+(𝑟_𝑏)²_𝑣+(𝑟_𝑏)²_𝑤. The position components of the perturbing body in the spacecraft topocentric frame are 𝑟_𝑏𝑢=𝑥_𝑏cos𝜗cos𝜑+𝑦_𝑏sin𝜗cos𝜑+𝑧_𝑏𝑟sin𝜑,_𝑏𝑣=−𝑥_𝑏sin𝜗+𝑦_𝑏𝑟cos𝜗,_𝑏𝑤=−𝑥_𝑏cos𝜗sin𝜑−𝑦_𝑏sin𝜗sin𝜑+𝑧_𝑏cos𝜑. (12) The perturbing acceleration is thus written as a function of time and state variables (viz., only 𝑟, 𝜗, and 𝜑, as gravity forces only depend on position). The lunisolar perturbation is the sum of the gravitational perturbations due to Moon (𝑏=𝑙) and Sun (𝑏=𝑙).

2.3. Solar Radiation Pressure

The photon pressure at distance 𝑅 from the Sun is 𝑝=𝐿_𝑆/4𝜋𝑅²𝑐_light, where 𝐿_𝑆 is the total power radiated by the Sun and 𝑐_light is the speed of light; the photon pressure at 𝑅^∗ = 1 AU is 𝑝^∗=4.55682⋅10⁻⁶N/m². Assuming reflectivity 𝜂=0.7 the acceleration on a spherical body of mass 𝑚 and cross-section 𝑆 is 𝐚_srp=1+𝜂𝑝^∗𝑅^∗𝑅²𝑆𝑚𝐑𝑅=Γ𝐑𝑚𝑅³, (13) whose components are 𝑎_srp𝑢=Γ𝑚𝑅³𝑟_𝑠𝑢,𝑎−𝑟_srp𝑣=Γ𝑚𝑅³𝑟_𝑠𝑣,𝑎_srp𝑤=Γ𝑚𝑅³𝑟_𝑠𝑤. (14) The effect of solar radiation pressure is therefore an acceleration in the Sun-spacecraft direction, inversely proportional to the squared distance of the two bodies. This acceleration and the solar gravity acceleration show the same dependence on distance and are parallel but with opposite directions; the similarity with the first term on the right-hand side of (11) allows one to treat them simultaneously. One should note that the perturbing acceleration in (13) depends also on the instantaneous mass; this fact introduces an additional term in the time derivative of the mass adjoint variable.

A conical shadow of the Earth is considered to determine the eclipses when (𝑟_𝑠)_𝑢<0 (Sun and spacecraft on opposite sides with respect to the Earth). The relevant quantities are sketched (not to scale) in Figure 2. The Earth determines a shadow cone with semiangle 𝛾_shadow=sin⁻¹(𝑟_𝐸/𝑟_𝑠), where 𝑟_𝐸 is the Earth radius. The spacecraft is instead on the surface of a cone (centered at the Sun with axis on the Earth-Sun connecting line) with semi-angle 𝛾=sin⁻¹(𝑟sin𝛿/𝑅), where the angle 𝛿 between the Earth centered vectors pointing to the Sun and spacecraft is evaluated as 𝛿=cos⁻¹[(𝑟_𝑠)_𝑢/𝑟_𝑠]. The spacecraft is in the Earth shadow when (𝑟_𝑠)_𝑢<0 and 𝛾<𝛾_shadow. Results show that, in the present case (𝑚=960kg, 𝑆=5.7m²), solar radiation pressure has a negligible influence on performance, causing variations of the final mass of a few grams.

Figure 2 

Schematic geometry of the Earth shadow.

3. Indirect Optimization Method and Optimal Controls

Optimal control theory (OCT) [29] is applied to the problem described above to determine the optimal solution. The state variables are collected in vector 𝐱=(𝑟𝜗𝜑𝑢𝑣𝑤𝑚)^𝑇. Adjoint variables 𝝀=(𝜆_𝑟𝜆_𝜗𝜆_𝜑𝜆_𝑢𝜆_𝑣𝜆_𝑤𝜆_𝑚)^𝑇 are introduced, and the Hamiltonian is defined as 𝐻=𝝀^𝑇𝐱=𝐻^′+𝚲^𝑇𝐓−𝜆_𝑚𝑇𝑐, (15) where 𝐻′ collects all the terms that do not contain the control (i.e., 𝐓in the present problem), and 𝚲=𝜆_𝑢𝜾+𝜆_𝑣𝚥+𝜆_𝑤𝐤 (16) is the adjoint vector to the velocity, which is named primer vector in the literature [30].

OCT provides the Euler-Lagrange equations for the adjoint variables d𝝀^𝑇d𝑡=−𝜕𝐻𝜕𝐱. (17) The perturbing terms depend on the state variables and therefore influence the derivatives of the adjoint variables. As far as the geopotential is concerned, differentiation with respect to 𝑟 and 𝜗 is again straightforward. The recursive scheme, which is outlined in Section 2.1, is again used to evaluate the derivatives with respect to declination 𝜑. Explicit expressions are derived for the gravitational perturbations of Moon and Sun (which depend on 𝑟, 𝜗 and 𝜑) and solar radiation pressure (which also depends on 𝑚). The equations of the derivatives, which are tedious but simple to obtain, are omitted for the sake of conciseness.

Pontryagin’s maximum principle (PMP), which states that the optimal control must maximize the Hamiltonian, is used to determine the optimal controls (thrust magnitude and direction). The thrust 𝐓must be parallel to the primer vector and the Hamiltonian, which is rewritten as 𝐻=𝐻^′Λ+𝑇𝑚−𝜆_𝑚𝑐, (18) is linear with respect to 𝑇; a bang-bang control arises and the thrust magnitude must be either maximum or minimum (i.e., zero) depending on the sign of the switching function 𝑆_𝐹=Λ/𝑚−𝜆_𝑚/𝑐 (Λ indicates the primer vector magnitude) 𝑇𝑇=^maxfor𝑆_𝐹>0,0for𝑆_𝐹<0. (19)

When indirect procedures are used, the thrust level is usually decided during integration on the basis of the sign of the switching function. It is instead convenient to split the trajectory into 𝑓 arcs that are joined at the switch points, that is, where the control (here the thrust magnitude) is discontinuous; in other problems a similar split is enforced where state variables are discontinuous and/or their values are subject to constraints. The time lengths of these arcs become 𝑓 additional unknown parameters of the present problem. This scheme avoids integration instabilities, assures high accuracy, and improves the convergence of the numerical procedure. The switching structure (i.e., number and order of thrust and coast arcs) is specified a priori; after the achievement of the numerical solution, the switching function history is checked and the switching structure is modified if PMP is violated (except at departure and arrival where, according to mission constraints, the engine cannot be used even though it would be beneficial).

OCT also provides the boundary conditions for optimality [3], which depend on performance index and boundary conditions on the state variables. Mayer formulation is adopted, and the performance index is written as 𝜙𝐱_(𝑗−1)+,𝐱_𝑗−,𝑡_(𝑗−1)+,𝑡_𝑗−𝑗=1,…,𝑓 (20) with subscripts 𝑗− and 𝑗+ indicating the values just before and after point 𝑗. Boundary conditions on the state variables are collected in the form 𝝍𝐱_(𝑗−1)+,𝐱_𝑗−,𝑡_(𝑗−1)+,𝑡_𝑗−=0,𝑗=1,…,𝑓. (21)

Optimality requires −𝝀^𝑇_𝑗−+𝜕𝜙𝜕𝐱_𝑗−+𝝁^𝑇𝜕𝝍𝜕𝐱_𝑗−=0,𝑗=1,…,𝑓, (22)𝝀^𝑇_𝑗++𝜕𝜙𝜕𝐱_𝑗++𝝁^𝑇𝜕𝝍𝜕𝐱_𝑗+=0,𝑗=0,…,𝑓−1, (23)𝐻_𝑗−+𝜕𝜙𝜕𝑡_𝑗−+𝝁^𝑇𝜕𝝍𝜕𝑡_𝑗−=0,𝑗=1,…,𝑓, (24)−𝐻_𝑗++𝜕𝜙𝜕𝑡_𝑗++𝝁^𝑇𝜕𝝍𝜕𝑡_𝑗+=0,𝑗=0,…,𝑓−1. (25) The constant Lagrange multipliers 𝝁are eliminated from (22)–(25); the resulting boundary conditions for optimality and the boundary conditions on the state variables, given by (21), are collected in a single vector in the form 𝝈𝐱_(𝑗−1)+,𝐱_𝑗−,𝝀_(𝑗−1)+,𝝀_𝑗−,𝑡_(𝑗−1)+,𝑡_𝑗−=0,𝑗=1,…,𝑓, (26) which, together with state and adjoint differential equations, defines a multipoint boundary value problem (MPBVP).

In the present problem, at the initial point (𝑗=0), 𝑡₀ and all state variables are assigned; at the final point (𝑗=𝑓) the apogee radius 𝑟_𝐴 and orbit semilatus rectum 𝑝 are given, and 𝑟_𝑓−𝑟_𝐴𝑢=0,_𝑓,𝑣=0_2𝑓+𝑤²_𝑓−𝜇𝑝𝑟²_𝐴=0 (27) are the conditions imposed on the state variables. The performance index to be maximized is the final mass, that is, 𝜙=𝑚_𝑓. The boundary conditions for optimality provide 𝜆_𝜗𝑓𝜆=0,_𝜑𝑓𝜆=0,_𝑣𝑓𝑤_𝑓−𝜆_𝑤𝑓𝑣_𝑓𝜆=0,_𝑚𝑓=1. (28)(29)(30)(31) The final time is free and (24) with 𝑗=𝑓 provides the transversality condition 𝐻_𝑓=0. (32)

Application of (24) and (25) at every switch point prescribes the Hamiltonian continuity; state and adjoint variables are continuous and, as a consequence, the switching function must be null at the switch points 𝑆_𝐹𝑗=0𝑗=1,…,𝑓−1. (33)

The numerical problem consists of 14 differential equations represented by (6) and (17). Initial values of the seven adjoint variables and time lengths of 𝑓 coast and burn arcs are unknown in the present problem; an equal number of boundary conditions, given by (27)–(33), completes the MPBVP. The problem is homogeneous in the adjoint variables, and (31) can be replaced by assigning the initial value 𝜆_𝑚0=1 in order to reduce the number of unknowns. The unknown parameters are collected in a vector 𝐩.

4. Numerical Solution

The solution of the problem outlined in the previous section is a difficult task; many aspects of the relevant numerical procedure deserve attention. Variables are normalized using the Earth equatorial radius, the corresponding circular velocity, and the spacecraft initial mass as reference values. Differential equations are integrated by using a variable-order variable-step integration scheme, based on the Adams-Moulton formulas, in order to carry out a fast and very accurate integration.

A single-shooting technique is adopted to solve the MPBVP, which is transformed into a series of initial value problems (IVPs) leading to convergence by means of Newton’s method [31], according to the following shooting algorithm: at the 𝑟th iteration (1) using the current tentative values 𝐩^𝑟 for the unknown parameters, the related IVP is solved numerically; (2) the errors 𝝈on the boundary conditions are evaluated; (3) the new tentative values 𝐩^𝑟+1 are evaluated according to Newton’s rule 𝐩^𝑟+1=𝐩^𝑟−𝜕𝝈𝜕𝐩⁻¹𝝈. (34) This process is repeated until all errors are smaller than an assigned value (here set equal to 10⁻⁷).

At any step, the sensitivity matrix 𝜕𝝈/𝜕𝐩is numerically evaluated by means of a first-order forward finite-difference scheme. Each unknown parameter is in turn perturbed by a small amount 𝛿𝑝_𝑖; the new IVP is solved, and the change of the errors on the boundary conditions 𝛿𝝈is evaluated. The choice of the perturbation step 𝛿𝑝_𝑖 is important to achieve a proper evaluation of the sensitivity matrix. A rule of thumb suggests 𝛿𝑝_𝑖 of the same order of the square root of the absolute tolerance used by the integrator (𝛿𝑝_𝑖=5⋅10⁻⁶ in this application).

The initial guess 𝐩¹ is very important for the convergence of the numerical procedure. The time at the switch points (i.e., where the engine is turned on or off) may exhibit significant changes (the orbital period is modified by the maneuvers that are performed). The corresponding right ascension instead can be easily guessed as thrust arcs are in correspondence of the orbit apsides, either perigee or apogee, and 𝜗 is preferred as the independent variable in place of time 𝑡. Moreover, 𝜗 is replaced by a nondimensional variable 𝜀 in order to “fix’’ the integration intervals; inside each arc 𝜀=𝑗−1+𝜗−𝜗_𝑗−1𝜗_𝑗−𝜗_𝑗−1, (35)𝜀 is zero at the initial point, assumes consecutive integer values at the switch points, and is equal to the total number of arcs 𝑓 at the final point. The 𝑓 right-ascension values 𝜗_𝑗 at the switch points and at the final point are the problem additional unknowns that replace the corresponding unknown times (𝜗₀=1.5𝜋 is an assigned initial condition).

The features of the present indirect approach, in particular, the peculiar treatment of the switch conditions, widen the convergence radius of the numerical procedure. On the other hand, the MPBVP solution is correct only if the tentative switching structure corresponds to that of the optimal solution. If the tentative solution lacks a burn arc, a suboptimal solution is usually found, with the switching function greater than zero during a coast arc. The request of adding a further thrust arc is clear, but, unfortunately, this sometimes implies the ensuing removal of another propelled arc. If the tentative structure presents one more burn than the optimal solution, (33) cannot be fulfilled at the two extremities of the superfluous arc, and the code would provide an unfeasible solution (e.g., a burn arc with negative length) or the MPBVP solution process would not converge at all. The suggestion that an arc should be removed is clear, but no replacement arc is necessary, and this case is easier to manage.

An experienced user is able to guess the unknown parameters and a suitable burn structure with a trial and error process, which is repeated until PMP is satisfied everywhere along the trajectory. This manual procedure is time-consuming, and an automated procedure is useful if the analysis has to cover a large number of different cases (e.g., a wide launch window). In the present problem the optimal switching structure is largely dependent on Moon and Sun positions and therefore on the departure date. A stepped procedure is devised to overcome the convergence difficulties related to the search for the optimal solution for an assigned departure date. The procedure exploits continuation techniques and is based on the observation that, in general, it is more difficult to introduce a new arc than to remove an existing one. It requires a tentative solution with the maximum number of thrust arcs that the problem could present. A suitable solution is obtained by assuming a simplified dynamical model that considers only the gravitational perturbation related to the second zonal harmonic; the optimal solution for this problem (termed J2 problem in the following) is independent of the departure date.

4.1. J2 Problem

The procedure starts using a dynamical model that considers only the gravitational perturbation related to the Earth oblateness, that is, 𝐶₂₀ is the only nonzero harmonic coefficient in (7). This problem is independent of the departure date.

According to the Keplerian model, the initial conditions of the spacecraft would permit ballistic attainment of the desired final apogee and the deployment would require only the perigee raising. The effect of the Earth oblateness on semimajor axis and eccentricity is null after a complete orbit; the perigee is unchanged, but the actual apogee of the initial orbit is lower. Actual differences between final and initial orbit in terms of perigee and apogee radii are 14650 km and 9172 km, respectively, when J2 is considered. The mission requires to raise both perigee and apogee, but the apogee maneuver is very expensive from a propulsive point of view, in comparison to the very small perigee maneuver. The optimal impulsive strategy [32] prescribes two burns: a perigee impulse, followed by an apogee impulse (PA is used to describe this burn sequence).

In the finite-thrust case it is convenient to split the impulsive maneuvers into multiple burn arcs centered at the apsides, in order to reduce the propulsive losses. The perigee burn, which should precede all the apogee burns, is very short and, if the number of revolutions is limited, the split of the longer apogee burn is preferable in terms of propellant consumption. The optimal burn structure consists of a single perigee burn (to raise the apogee) followed by burns at every apogee passage. When 4.5 revolutions are permitted, only three apogee burns can be performed (engine cannot be used in the proximity of departure and arrival apsides): the optimal burn sequence is therefore PAAA.

After this suitable switching structure is assumed, one can easily guess the right ascension at the engine switches. The physical meaning of the primer vector suggests that the adjoint vector to velocity should be parallel to the thrust direction, which is essentially horizontal and in the orbital plane for the initial burn. Latitude has a minimal influence, and longitude has no influence at all; the corresponding adjoint variables are therefore (roughly) zero. From a practical point of view, the magnitude of the primer vector Λ and the adjoint variable to radius 𝜆_𝑟 are the only parameters difficult to estimate; convergence is easily obtained.

4.2. Stepped Procedure

When the complete dynamical model is considered, the optimal trajectories (and the associated control laws) are dependent on time and their burn structure becomes more difficult to assess. Moon and Sun with their attraction are capable of varying significantly the perigee radius; the Moon, in particular, can raise/lower the perigee even by 200 km in a single ballistic revolution, depending on its position. Thus, the optimal length of the apogee burns is dictated by the need not only of containing the propulsive losses but also of modifying the orbital periods; the time of passage at each consecutive apogee is anticipated or delayed, in order to find the perturbing bodies in more favorable or less unfavorable positions.

As a consequence, the optimal solution of the complete problem does not exhibit a uniform split of the three apogee burns. The optimal split depends on the departure date: in many cases one burn (occasionally two burns) may vanish. The switching structure of the mission departing on an assigned date cannot be guessed a priori. A suitable procedure that autonomously fixes the optimal switching structure is here outlined. Starting from the solution of the J2 problem, which is independent of the departure date and has engine firings during all apogee passages, a fraction 𝑃_𝑓 of the remaining perturbations is introduced, and the solution of the corresponding MPBVP is searched for; the fraction 𝑃_𝑓 is progressively increased in a discrete number of steps (𝑘_𝑃) until it reaches unity.

At each step, 𝑃_𝑓 is not imposed as an assigned constant but is considered as an unknown parameter, with an additional constraint that enforces it to the desired value. The iterative procedure progressively moves 𝑃_𝑓 towards its target value, while the other unknowns are simultaneously adjusted to satisfy all the boundary conditions. In the first 𝑘_𝑅 iterations a relaxation of the Newton scheme, which uses only a fraction of the correction in (34), is used to avoid that abrupt changes of the unknown parameters prevent convergence.

Convergence is usually reached if the step Δ𝑃_𝑓 is chosen small enough. Solution is not achieved (or an unfeasible solution is found) only when the optimal switching structure has one burn less than the previous step. In these cases, the optimal solution, which was obtained with the previous level of 𝑃_𝑓, suggests which thrust arc has to be removed. In particular, the maximum values of the switching function during every burn arc are compared, and the burn arc with the smaller value is removed. The right ascensions at the extremes of the arc to be removed are imposed to coincide with the closer apsis (thus imposing null length of the burn arc); these new boundary conditions replace (33) at the same extremes. The last converged solution is still used as initial guess; at the first iteration the only unsatisfied boundary conditions concern perturbation fraction and right ascensions at the extremes of the arc to be removed. An analysis of the achieved solution in the light of PMP is performed to confirm that the new switching structure is optimal.

The suitable values for 𝑘_𝑃 and 𝑘_𝑅 are problem dependent; in general, one can expect to reduce issues related to convergence when these two values are increased, but this would also increase the overall computation time. For the problem under investigation a good balance is given by choosing 𝑘_𝑃=7 (with 𝑃_𝑓=0.05,0.1,0.2,0.4,0.6,0.8,1) and 𝑘_𝑅=4 (1/5, 1/4, 1/3, and 1/2, in this order, are the fractions used). The computational time to obtain a full-perturbation 4.5-revolution solution is about 30 seconds on an Intel i7 CPU at 2.67 GHz.

The procedure, which is fully automated, has been tested for 366 departure dates between December 1, 2015, and November 30, 2016, with one-day steps. An analysis of the results has shown that the optimal solution according to PMP is achieved almost in every case. In a single instance the procedure removed a very tiny arc that is present in the optimal solution. This error was due to insufficient values for 𝑘_𝑃 and 𝑘_𝑅; convergence was immediately obtained by adding 4 additional steps and increasing 𝑘_𝑅 to 9. Moreover, some of the solutions without the first apogee burn (P0AA) are only suboptimal. The corresponding optimal solutions present two ballistic revolutions, followed by a perigee and two apogee burns (0PAA); however, final masses differ only by a few grams. The correct solution can be achieved manually; an autonomous convergence procedure was not developed due to the scarce number and interest of these low-performance opportunities.

The applicability of this procedure is more general and may concern any other parameter whose progressive change reduces the number of burn arcs. For instance, an analogous procedure could be adopted to deal with changes of the thrust magnitude; in particular, when the thrust is low enough, the engine is on at every apsis passage, but burn arcs may vanish when the thrust magnitude is increased: the analysis of the switching function suggests which arc has to be removed. It was verified that an 8 N solution for the deployment here considered can be obtained from the 1 N solution departing on the same date, by progressively increasing the thrust magnitude (e.g., in 1 N steps).

5. Numerical Results

The deployment of the spacecraft into the target HEO was studied in details using the indirect procedure previously described. Perturbations (except J2) were initially neglected, in order to assess the general characteristics of the optimal strategy. Two thrust levels were considered with the aim of better understanding how thrust acceleration and number of revolutions influence the propellant consumption. No burns are permitted in the first revolution and propulsion cannot be employed at the final apogee insertion. The maneuver requires at least 2.5 revolutions around the Earth; a maximum of 4.5 revolutions was assumed.

5.1. J2 Solution

Figure 3 shows that the final mass increases with the number of revolutions and burns; this effect is more important when the thrust acceleration is low. In fact, the number of burns has no influence for impulsive thrust, as impulses can be applied exactly at apogee and perigee passages, without velocity losses; these instead arise for the finite-thrust maneuvers as thrust is exploited far from the apsides. A larger number of revolutions permit the split of the propulsive effort between shorter arcs centered at the apsides, with lower losses. One should also remember that a larger number of burns are preferable as they allow an easier correction of the errors that may occur during the maneuvers.

Figure 3 

Final mass for different numbers of revolutions (J2 only).

Table 2 provides some details of these maneuvers. The different time length of the thrust arcs during the same transfer is related to the mass consumption; as the spacecraft becomes lighter, thrust arcs are shortened in order to have an almost uniform split of the total velocity change Δ𝑉 between the burns. Note that the velocity is larger at each later apogee and the spacecraft sweeps a longer arc even though the corresponding time length is shorter. This effect is again more evident for the lowest thrust level. One should also note that the time length of the perigee burn increases with the number of revolutions. This happens because, during the long apogee arc (A1) of a fast transfer, the apogee height is slightly increased, reducing the requirement for the perigee burn, as shown in Figure 4, which presents the evolution of perigee and apogee for the 2.5-revolution transfer with 𝑇=8 N. The apogee change during the apogee burns (A1, A2, and A3) when 4.5 revolutions are performed is instead negligible, as shown in Figure 5. The 𝑟_𝑎 and 𝑟_𝑝 variations, which occur during the coast arcs, are due to the J2 effect.

Table 2 

Time length and angular length of burn arcs (J2 only).

Figure 4 

Perigee and apogee evolution for the 2.5-revolution 8 N transfer (J2 only).

Figure 5 

Perigee and apogee evolution for the 4.5-revolution 8 N transfer (J2 only).

The shortest transfer with 𝑇=1 N is the only maneuver with two perigee burns as the time length of the apogee firing would otherwise not be sufficient to obtain the required perigee rise. For this reason, two perigee burns are necessary: the first one (P1) increases the apogee above the required value (see Figure 6) to have a longer available time for an efficient perigee rise during the apogee burn; the second perigee burn (P2) eventually reduces the apogee to the required value.

Figure 6 

Perigee and apogee evolution for 2.5-revolution 1 N transfer (J2 only).

5.2. Lunisolar Perturbation

Perturbations other than J2, mainly solar and lunar gravity, change the performance of the deployment maneuver and the optimal switching structure. The final mass for the 4.5-revolution transfer with 𝑇=8N is shown in Figure 7 as a function of the departure date in the 1-year launch window starting December 1, 2015 (MJD 57357), which is studied in detail. Figure 8 shows the final mass in the first three months of the launch window and indicates the burn structure of each optimal transfer.

Figure 7 

Final mass for 4.5-revolution 8 N transfers.

Figure 8 

Switching structure for 4.5-revolution 8 N transfers.

A simplified analysis permits an estimation of the most significant effect of the third body gravitation, which affects the spacecraft perigee height, as a function of the perturbing body position. Since the orbit is highly elliptic, the spacecraft spends most time at the apogee (𝜈=180deg, 𝜗=90deg) where, in addition, the perturbing acceleration is larger due to the greater Earth-spacecraft distance. Therefore, only this point is considered in the following analysis. An increase of the perigee altitude is the mission main requirement, and the acceleration component parallel to the apogee velocity, that is, the tangential component 𝑎_𝑡, is the main cause of the investigated effect. If coplanar orbits are assumed, one easily determines (see Figure 1) 𝑎_𝑡𝜇=−^𝑏𝑟²_𝑏𝑟_𝑏𝑅³−1cos𝜗_𝑏 (36) with the spacecraft distance from the perturbing body expressed as 𝑅²=𝑟²_𝑏+𝑟²−2𝑟_𝑏𝑟sin𝜗_𝑏. (37)

When the Sun is the perturbing body, 𝑟≪𝑟_𝑏 and only first-order terms are retained to obtain (𝑟_𝑏/𝑅)³≈1+3(𝑟/𝑟_𝑏)sin𝜗_𝑏 and 𝑎_𝑡𝜇≈−^𝑏𝑟²_𝑏32sin2𝜗_𝑏 (38) with maximum positive values at 𝜗_𝑏=135 and 315 deg (the most favorable positions of the Sun) and maximum negative values at 𝜗_𝑏=45 and 225 deg (the most unfavorable positions). If short-period oscillations due to the lunar gravity are ignored (two cycles are found in each sidereal period of the moon), the final mass in Figure 7 closely follows this trend, with two sinusoids in the 1-year launch window. The peak values occur for departure on January 24, 2016, (MJD 57411, right ascensions of the Sun at departure and arrival are 307 and 323 deg) and July 27, 2016 (MJD 57596, 127 and 143 deg).

When the Moon is considered, the spacecraft distance from the Earth becomes comparable to the Earth-Moon distance (𝑟/𝑟_𝑏≈0.5) and the previous simplification does not hold. The symmetry of the result with respect to the 𝑥 axis is broken and the effects of the third-body perturbation are enhanced when spacecraft apogees and Moon are on the same side with respect to the Earth, that is, when sin𝜗_𝑏>0. The maximum benefit occurs when 𝜗_𝑏≈115deg (with a less pronounced beneficial effect at 𝜗_𝑏≈330deg), whereas the largest negative effect is at 𝜗_𝑏≈65deg (with a less remarkable effect at 𝜗_𝑏≈210deg). Figure 9 shows the performance in terms of final mass for the 2.5- and 4.5-revolution transfers, when the Moon is held fixed during the whole maneuver and the solar attraction is neglected. The longest missions exploit significant lunar assists during five apogee passages, two more than in the shortest missions: variations of the final mass with respect to the average value show roughly the same 5/3 ratio.

Figure 9 

Final mass as a function of the Moon fixed position for 2.5- and 4.5-revolution 8 N transfers (Sun neglected).

The moon influence becomes more complex when its actual motion is considered. The final mass is shown in Figure 10 as a function of the Moon position at the beginning of the maneuver; the solar gravitational attraction is again neglected. In this case, the final mass presents smaller variations than in the fixed-Moon case (𝑦-axis scale is different from Figure 9). During each spacecraft revolution the Moon moves about 40–55 degrees (changes are due to the Moon eccentricity and to the increase of the spacecraft orbital period caused by apogee burns) and will approximately occupy an unfavorable position two revolutions after a favorable configuration and vice versa.

Figure 10 

Final mass as a function of the Moon initial position for 2.5- and 4.5-revolution 8 N transfers (Sun neglected).

The Moon is at the most favorable position (𝜗_𝑏=115deg) at the first apogee passage of the best 2.5-revolution transfer, offsetting the penalty at the last one (205 deg). On the contrary, the worst performance occurs when the three apogee passages find the Moon at about 20, 70, and 125 degrees, with the most unfavorable configuration at the second passage.

The Moon moves about 210 degrees during the best 4.5-revolution transfer and is in a favorable position at the first (115 deg) and last (about 305 deg) apogee passages (with a single unfavorable position, at about 205 degrees, during the third passage). However, the spacecraft adjusts the burn lengths and varies the orbital period during each revolution to put forward or push back the passages in order to enhance/reduce the effects of favorable/unfavorable geometrical configurations. On average, the final mass is larger in comparison to the shortest transfers, but the maximum achievable mass is lower.

The switching structure of the optimal 4.5-revolution missions (see Figure 8) changes with a clear regularity according to the departure date. The best missions require the removal of the last burn arc (PAA0), whereas the removal of the first arc (P0AA) is beneficial in the worst cases. Optimization prescribes the removal of the last two arcs in four cases, and the switching structure becomes PA00. The same optimal structures repeat at roughly 14-day intervals, corresponding to a half revolution of the Moon around the Earth. For an assigned departure date, the thrust strategy has almost no capability of phasing the initial apogee passage with the Moon position. However, a longer thrust arc at the first apogee increases the total time of flight. On the contrary, when the first apogee burn vanishes, the following orbital periods are shorter and the whole mission is faster. The trip time may differ more than 12 hours (about 6 degree in angular position of the Moon). Figure 11 shows that the mission departing on December 12 delays the last apogee passage to find the Moon in a more favorable position. On the contrary, but with a similar aim, the mission starting on December 5 anticipates the fourth apogee passage. The Moon complex influence on the spacecraft trajectories suggests that further analyses could provide interesting hints for all the missions that exploit lunar resonance.

Figure 11 

Moon angular position during 4.5-revolution transfers.

5.3. Example Case Analysis

The deployment with departure on December 1, 2015, is analyzed in detail to describe the continuation procedure, which progressively increases the perturbation magnitude taken into account. Figure 12 shows an enlargement of the switching function history as a function of the nondimensional variable 𝜀, for different levels of perturbation fraction 𝑃_𝑓; this representation is chosen for the sake of clarity, even though it hides the actual arc lengths (provided in Table 3). Thrust would be required at departure and arrival where PMP is not satisfied (𝑆_𝑓>0) but is not allowed for operational reasons. The odd intervals correspond to coast arcs, whereas the second arc is the perigee burn and the other even intervals correspond to three apogee burns.

Table 3 

Characteristics of the 4.5-revolution transfer with departure on December 1, 2015.

Figure 12 

Switching function history for 4.5-revolution transfer departing on December 1, 2015.

At the beginning of the stepped procedure, when no perturbation except J2 is considered (𝑃_𝑓=0), the switching function is almost the same during every apogee burn. For this particular departure date, when the lunisolar perturbation is increased, the first apogee burn becomes more convenient, and in this arc the switching function peak grows, while the burn length increases; on the other hand, the last thrust arc is shorter, and there the switching function peak decreases. The switching structure remains the same until 𝑃_𝑓=0.6 but provides an unfeasible solution for 𝑃_𝑓=0.8 and must be changed. The lowest peak of the switching function for 𝑃_𝑓=0.6 occurs in the third apogee burn, which is therefore removed. Convergence assuming the new PAA0 structure is obtained, and PMP confirms the optimality of this solution. Arc time lengths and relevant masses for the December 1 departure are summarized in Table 3.

Figure 13, which refers to the deployment with departure on December 1, 2015 (8 N thrust, 4.5 revolutions), shows the angle 𝛽 between thrust and orbit plane (positive values tend to increase inclination) and the angle 𝛼 between the spacecraft velocity and the thrust projection onto the orbit plane (positive values correspond to thrust towards the Earth). When all perturbations are considered, a small out-of-plane thrust component is introduced during each burn in order to slightly change the orbit plane and improve the beneficial effect of the luni-solar gravitational perturbation. In the orbit plane, the thrust is directed inward with respect to the Earth when the spacecraft is moving outward (0<𝜈<180deg, i.e., after perigee and before apogee) and vice-versa, to reduce the orbit eccentricity in agreement with the variational equations for orbital parameters.

Figure 13 

Thrust angles of the 4.5-revolution transfer with departure on December 1, 2015.

6. Conclusions

The maneuver for the deployment of a satellite into a highly elliptic orbit has been analyzed by means of an indirect optimization method, which handles perturbations from Earth asphericity, lunar and solar gravity, and solar radiation pressure. The optimal thrust strategy depends greatly on the departure date, and is difficult to assess “a priori,’’ especially for multirevolution missions. A continuation scheme based on the gradual introduction of the relevant perturbations is used to converge to the optimal solutions starting from a single initial solution, obtained only considering J2 effect. The switching function is checked during the procedure to maintain the optimal burn structure. The algorithm proves to be fast and reliable and convergence is almost always obtained automatically, without any user’s action. Simple adjustments of the procedure parameters may be necessary when, in rare cases, convergence to the optimal solution is not directly obtained.

The availability of a powerful tool permitted a detailed analysis of this kind of deployment. The influence of thrust level and number of revolutions on the transfer performance, that is, mass delivered to the final orbit, has been investigated. The benefit of an increased number of burns and the performance drop that is experienced when a fast transfer is combined with low thrust level have been quantified. As far as the lunisolar perturbation is concerned, the position of the perturbing bodies in correspondence of the spacecraft passages at apogee has a major influence on performance. Due to its slow apparent motion, the Sun position that gives the maximum benefit can be exploited by properly choosing the launch date, without significant changes of the thrust strategy, that is, of the switching structure. On the contrary, the Moon moves rapidly and occupies different positions at the consecutive spacecraft apogee passages; in the search for the optimal maneuver, the time lengths of coast and thrust arcs change (in limit cases, some burn arcs may disappear and the optimal switching structure is modified) in order to find the Moon in the most favorable positions. Results show that the proper selection of the departure date and the use of an optimal strategy, which exploits luni-solar perturbations, can save several kilograms of propellant.

Acknowledgment

This work was supported by the Centre National d’ Etudes Spatiales (CNES; Contract no. 93333/00).