Optimal Space Trajectories with Multiple Coast Arcs Using Modified Equinoctial Elements

The detection of optimal trajectories with multiple coast arcs represents a significant and challenging problem of practical relevance in space mission analysis. Two such types of optimal paths are analyzed in this study: (a) minimum-time low-thrust trajectories with eclipse intervals and (b) minimum-fuel finite-thrust paths. Modified equinoctial elements are used to describe the orbit dynamics. Problem (a) is formulated as a multiple-arc optimization problem, and additional, specific multipoint necessary conditions for optimality are derived. These yield the jump conditions for the costate variables at the transitions from light to shadow (and vice versa). A sequential solution methodology capable of enforcing all the multipoint conditions is proposed and successfully applied in an illustrative numerical example. Unlike several preceding researches, no regularization or averaging is required to make tractable and solve the problem. Moreover, this work revisits problem (b), formulated as a single-arc optimization problem, while emphasizing the substantial analytical differences between minimum-fuel paths and problem (a). This study also proves the existence and provides the derivation of the closed-form expressions for the costate variables (associated with equinoctial elements) along optimal coast arcs.


Introduction
Spacecraft trajectory optimization is concerned with the determination of the optimal direction and magnitude of the propulsive thrust that drive a space vehicle toward some specified conditions, while minimizing either propellant consumption or time of flight. Optimal space trajectories have been extensively investigated from both the analytical and the numerical perspective, using a variety of approaches. The impulsive thrust assumption [45] represents an excellent approximation for spacecraft that employ chemical propulsion for short durations. In this context, orbit transfer optimization, aimed at minimizing propellant consumption, reduces to a nonlinear programming problem. However, in the presence of moderate or low thrust levels, the impulsive approximation loses accuracy, and the general properties of optimal finite-thrust trajectories can no longer be inferred from an impulsive solution [34,49]. In these cases, the optimal path of interest must be found as the solution of a continuous-time optimal control problem.
As a valuable alternative to chemical propulsion, in recent years low-thrust electric propulsion [48] attracted an increasing interest by the scientific community, and already found application in a variety of mission scenarios, e.g. the NASA Deep Space 1 and the ESA Smart-1 missions [43,44]. Thanks to high values of the specific impulse, low-thrust propulsion allows substantial propellant savings, at the price of increasing-even considerably-the time of flight. Pioneering studies on low-thrust trajectories are due to Edelbaum [13], who apparently was the first scientist to point out the advantages of using low-thrust in space missions. Most recently, extensive research on the same subject was carried out by Petropoulos [35,36], Betts [5][6][7], Ross [46], and Kechichian [23][24][25][26][27], to name a few. Low-thrust trajectory optimization problems are solved through the use of direct, indirect, or heuristic approaches sometimes combined in hybrid forms [37]. With this regard, Betts [8], Rao [42], and Conway [11] offer excellent overviews of the available methods in spacecraft trajectory optimization. A major drawback of low-thrust electric propulsion resides in the demanding requirement of electrical power to operate. As a result, in operational scenarios low-thrust propulsion is switched off when the satellite is eclipsed. It is quite obvious that this circumstance implies a complete reconsideration of the underlying optimization problem. In general, this is formulated as a minimum-time problem, for which the available control depends on the state. A limited number of studies are devoted to low-thrust trajectory optimization with eclipsing. The cylindrical shadow model is assumed in several works as a very accurate approximation, due to the distance that separates the Sun from the Earth [15,17]. Under some simplifying assumption, Kechichian [26,27] derived analytical expressions for the variation of the orbit elements due to low thrust, with the inclusion of the eclipse effect on the availability of electrical power. Some researchers employed orbital averaging, in conjunction with either direct techniques [16,28], or the sequential gradient-restoration algorithm [2], or a shooting method [15]. Recently, Kluever [29] provided an algorithm, based on Edelbaum's analytic solution, devoted to computing the transfer time for low-thrust maneuvers with eclipsing, while Betts [7] proposed a direct optimization method that includes a multiple-phase First, virtually all types of trajectories can be described using modified equinoctial elements, unlike what occurs if the classical orbit elements are employed. Second, 5 out of 6 equinoctial elements remain constant (while the sixth is integrable) along coast arcs, in the presence of a single attracting body. In more complex mission scenarios, the space vehicle may be subject to perturbing accelerations, other than the action of the dominating attracting body, and in this case these 5 elements exhibit slow time variations. Third, in the numerical solution of low-thrust path optimization problems (with no eclipse intervals), the use of equinoctial elements was proven to mitigate the hypersensitivity issues encountered with spherical coordinates [38,47]. In fact, for multi-revolution orbit transfers, they exhibit superior convergence properties if compared to spherical coordinates, and are also amenable to homotopy methods [48]. With reference to problem (a), the present study has the primary objectives of (1) formulating the problem as a multiple-arc optimization problem, (2) deriving the complete set of necessary conditions for optimality, (3) proving that all the multipoint conditions referred to intermediate times can be solved sequentially in the numerical integration process, and (4) employing the latter conditions in a numerical example, for illustrative purposes. No averaging or regularization is introduced to make the problem tractable. Then, this research revisits problem (b) using equinoctial elements, with the objective of deriving the necessary conditions associated with minimum-fuel paths, while emphasizing the substantial analytical differences from problem (a). As a further contribution, the optimal adjoint equations along coast arcs are analyzed, with the intent of ascertaining the existence of closed forms for the costate variables associated with modified equinoctial elements.

Orbit Dynamics
This research considers a space vehicle that orbits a single celestial body, in the dynamical framework of the restricted problem of two bodies. The spacecraft of interest is modeled as a point mass.
In general, orbital motion can be described using either Cartesian coordinates, spherical variables, or osculating orbit elements, i.e. semimajor axis a, eccentricity e, inclination i, right ascension of the ascending node (RAAN) Ω , argument of periapse , and true anomaly f [41]. However, the Gauss equations [41], which govern the time evolution of the orbit elements, become singular in the presence of a circular or equatorial orbit (and also when an elliptic orbit transitions to a hyperbola). For these reasons, the modified equinoctial elements [9] were introduced, and are selected in this work as the variables that identify the dynamical state of the space vehicle. These elements are defined as [3,9] It is straightforward to recognize that x 1 represents the orbit semilatus rectum. Unlike the classical orbit elements, the modified equinoctial elements are never singular, (1) x 2 = e cos (Ω + ) x 3 = e sin (Ω + ) with the only exception of i = (condition that is unlikely to encounter, because equatorial retrograde orbits are rather impractical). If ∶= 1 + x 2 cos x 6 + x 3 sin x 6 , the instantaneous radius is r = x 1 ∕ . The classical orbit elements can be retrieved by inverting Eq. (1). The spacecraft position can be written in terms of a, e, i, Ω , , and f [41] or can be computed directly from the equinoctial elements [19].

Equations of Motion
The dynamical evolution of the modified equinoctial elements is governed by the respective Gauss equations [3], where is the gravitational parameter of the attracting body. The terms a r , a , a h are the components of the non-Keplerian acceleration a in the local vertical local horizontal (LVLH) frame aligned with r,̂,ĥ , where unit vector r is directed toward the instantaneous position vector r (taken from the center of the attracting body), whereas ĥ is aligned with the spacecraft angular momentum. The modified equinoctial elements allow identifying the instantaneous position and velocity of the spacecraft. This is controlled using the thrust supplied by the propulsion system. Let T max and m 0 represent the maximum available thrust magnitude and the initial mass of the space vehicle. If x 7 denotes the mass ratio and T the thrust magnitude, for x 7 the following equation can be obtained: where c represents the (constant) effective exhaust velocity of the propulsion system, whereas m is the instantaneous mass. The magnitude of the instantaneous thrust acceleration is a T = u T m 0 m = u T ∕x 7 and is constrained to the interval In the end, the spacecraft dynamics is described using the state vector x and the control vector u defined as In light of Eq. (10), the state Eqs. (2)-(8) can be written in compact form as

Spacecraft Eclipse Condition
Some types of low-thrust propulsion systems [48] require onboard electrical power in order to operate. In most cases, this can be provided only when the space vehicle is illuminated. This implies that propulsion must be considered unavailable when the spacecraft is eclipsed. This subsection focuses on the eclipse condition for space vehicles that orbit the Earth.
As a preliminary step, the spacecraft position in a proper inertial frame is derived. The Earth-centered inertial frame (ECI) has origin at the Earth center and axes aligned with the right-hand sequence of unit vectors ĉ 1 ,ĉ 2 ,ĉ 3 . Vector ĉ 1 identifies the vernal axis, which corresponds to the intersection of the Earth equatorial plane with the ecliptic plane, whereas ĉ 3 points toward the Earth rotation axis. In the ECIframe, the spacecraft Cartesian coordinates can be written in terms of the classical orbit elements as shown in Ref. 41. Then, relatively long trigonometric steps (omitted for the sake of brevity), based on the definitions (1), lead to the following expressions for the three coordinates X, Y, Z: where c ∶= cos and s ∶= sin ( denotes a generic angle).
Under the very reasonable (approximating) assumption that the Earth describes a circular orbit about the Sun, the position of the latter in the ECI-frame is identified by the three coordinates X S , Y S , and Z S , where r S equals 1 AU, (= 23.4 deg) is the ecliptic obliquity, whereas S identifies the instantaneous position of the Sun. If S0 denotes the value of S at the initial time t 0 (set to 0) and S = 2 T S ; T S = 1 year is the (constant) angular rate of the Sun in its motion relative to Earth, then S = S t + S0 .
Once the spacecraft and Sun positions are specified in the ECI-frame, the condition for eclipsing is [12] where R E is the Earth radius. Because r S ≫ R E , 2 ≃ ∕2 , and the eclipse condition (14) becomes Insertion of Eqs. (12) and (13)

into (15) yields
Usefulness of the previous relation will be apparent in the next section.
It is worth remarking that some trajectory arcs can be traveled in penumbral or antumbral conditions, related to the relative positions of Earth, Moon, and Sun. However, in the great majority of the mission scenarios these arcs are so short that propulsion can remain ignited. For this reason, penumbral and antumbral conditions are neglected in this study.

Minimum-Time Trajectories with Eclipse Intervals
Low-thrust propulsion systems usually require considerable electrical power to operate. For this reason, they are ignited only when the satellite is illuminated. This section is devoted to investigating the minimum-time transfer paths from an initial to a final Earth orbit, assuming availability of the thrust only in the intervals when the space vehicle is illuminated.

Statement of the Problem
The spacecraft of interest is governed by the state equations (11) and is subject to some (problem-dependent) boundary conditions of the form where denotes a generic column vector, which depends on the final time t f and on the initial and final state 0 and f . These conditions usually include the relations that define the initial and final orbits. As an example, the two terminal orbits can be defined by the respective values of the first five equinoctial elements, and in this case includes 10 components of the form where the symbols with bar denote specified values, whereas subscripts 0 and f respectively refer to the initial and final value of the corresponding variable. The initial time is assumed specified and is set to 0. The objective function J to minimize is the time of flight, therefore For the orbit transfer problem at hand, the thrust is available only when the spacecraft is illuminated, i.e. when inequality (16) is violated. It is quite obvious that the entire trajectory is composed of an unspecified number N of eclipse arcs and light arcs. Therefore, this becomes a multiple-arc trajectory optimization problem. In each arc j, the state equations are Thrust is switched off along the eclipse arcs ( u (j) (16) is fulfilled. This means that the optimal control law must be determined only in the light arcs. Unlike the control vector u, the state x is continuous across two adjacent arcs, i.e.
For each variable the subscripts ini and fin denote the initial and final value in the arc with index reported in the superscript. Furthermore, the transition between two adjacent arcs (either an eclipse arc followed by a light arc or a light arc followed by an eclipse arc) corresponds to fulfillment of the equality associated with Eq. (16), rewritten in compact form as Journal of Optimization Theory and Applications (2021) 191:545-574 The symbol (j) fin collects all the state components that appear in the left-hand side of Eq. (16), evaluated at the end of arc j, which occurs at time t j . It is worth noticing that (j) depends explicitly on t j through the variable S , evaluated at t j i.e. S,j = S t j + S0 . In the end, the multiple-arc problem consists in finding the optimal control u that minimizes the objective function (18), while holding the state equations (19) and the multipoint conditions (17), (20), and (21). Multiple-arc optimal control problems were also investigated by the author in Ref. [39].

Necessary Conditions for Optimality
In order to derive the necessary conditions for optimality, the extended objective function J is defined, where , j , and j are time-independent adjoint variables conjugate to the multipoint conditions (17), (20), and (21), whereas (j) is the time-varying costate vector associated with the state equations (19) and t N ≡ t f . N Hamiltonian functions H (j) (j = 1, … , N) and a function Φ are introduced, and the extended objective function J is rewritten as The first differential dJ can be obtained by using the steps illustrated in Appendix, and is where subscripts 0 and f refer to the initial and final time, respectively, whereas the symbol denotes the variation, i.e. the time-fixed differential, using the notation of Ref. 22 (cf. also Appendix). The first differential must vanish at an extremal [22], for arbitrary values of ,…,N , and this implies that the following necessary conditions for optimality must hold: The last condition, which implies stationarity of H (j) with respect to (j) , can be replaced by the more general Pontryagin minimum principle [32], (23), whereas the subscript * denotes the optimal value of the corresponding variable. Equation (33) is the adjoint equation for the costate variables, accompanied by the related boundary conditions (30) and (31) and by the multipoint conditions (27) and (28). Because appears in Eqs. (30) and (31), their explicit form is problem-dependent, unlike what occurs for the adjoint equations (33). Moreover, Eqs. (29) and (32) hold for the Hamiltonian functions, evaluated at times t 1 , … , t N−1 , t N ≡ t f . While Eqs. (30)- (35) are formally identical to the necessary conditions that hold for ordinary (single-arc) optimization problems, Eqs. (27)-(29) represent additional relations, termed multipoint necessary conditions henceforth.
(2)- (9) and (23), H can be rewritten as where r , , h , and (whose expressions are not reported for the sake of brevity) The adjoint variable (j) 7 , which is the seventh component of (j) , deserves special attention. Due to Eqs. (33) and (36), the costate equation for (j) 7 is Moreover, because (j) is independent of x (j) 7 (cf. Eq. (16)), combination of Eqs. (27) and (28) yields Because the final mass is unspecified, the boundary conditions (17) are independent of x 7,f . As a result, Eq. (31) yields While in the eclipse arc the spacecraft is subject to no control due to unavailability of low-thrust, in light arcs the minimum principle (35) allows expressing the optimal control in terms of the state and costate variables. With reference to Eq. (36), the first three terms in square parentheses can be regarded as a dot product. Thus, since u T ∕x 7 ≥ 0 , the thrust angles that minimize H (j) are given by (35) (j) * = arg min (j) H (j) (j) Using these expressions for (j) and (j) , Eqs. (36) and (37) become The latter relation, in conjunction with the final condition (39) and the continuity condition (38), implies that (j) 7 is nonnegative at all times. Moreover, because  1, … , N) , , k , and k (k = 1, … , N − 1).

Sequential Solution of the Multipoint Necessary Conditions for Optimality
The formulation of the orbit transfer with eclipse intervals as a multiple-arc optimization problem leads to establishing an extended set of necessary conditions for optimality. In particular, the multipoint necessary conditions for optimality (27)- (29), together with the multipoint conditions (20) and (21), represent a considerable number of additional relations to satisfy in the numerical solution process. However, Eqs. (27)-(29) can be combined, for the purpose of developing an advantageous methodology to enforce them.
As a first step, insertion of Eq. (27) into (28) yields (41) (j) fin , because the state is continuous across adjacent arcs. Two cases can occur, with regard to Eq. where I denotes the identity matrix, with dimension appropriate to the context. Equation (47) can be regarded as a linear system with (j+1) ini as the unknown. Therefore, (j+1) ini can be obtained using Eq. (47), if all the other quantities (including (j) fin ) are known.  ini , (j) fin . In fact, as the eclipse arc duration reduces or as u (max) T tends to zero, (j+1) ini must tend to (j) fin , to retrieve continuity of the adjoint variables, which holds for single-arc problems with no discontinuity in the available thrust magnitude.
It is apparent that Eqs. (47)-(48) represent the matching conditions for the adjoint vector at times t j (j = 1, … , N − 1) . In the scientific literature, Cerf [10] provided the relations that identify the jump conditions for the adjoint variables associated with Cartesian coordinates, by employing the Weierstrass-Erdmann corner conditions. Equations (47)-(48) represent the jump relations for the adjoint variables associated with equinoctial elements. They are derived making reference to a transition condition that has the general form reported in Eq. (21), and using the general principles of variational calculus (cf. also Sect. 3.2 and Appendix), in conjunction with the Pontryagin minimum principle. It is straightforward to recognize that Eq. (38) is the seventh scalar equation associated with either Eq. (47) or Eq. (48). Based on these relations and the remaining necessary conditions for optimality, the numerical solution process can include the following steps: Step 1. Identify the known initial values of the state and costate variables using Eqs. (17) and (30).
Step 2. Derive the (possible) relations between the initial values of the state and costate variables (i.e., the components of ) by eliminating some components of from Eq. (30); this leads to identifying the minimal set of unknown initial values for the components of 0 and 0 .
Step 3. Select the unknown initial values for the state and costate components belonging to the minimal set identified at Step 2; calculate the remaining initial values of the state and costate components using the relations found at Step 2.
Step 4. Select the final time t f .
Step 5. Until the current time t ≤ t f , iterate the following sub-steps: identify the type of arc (either a light arc or an eclipse arc); integrate numerically the state and costate equations (19) and (33), setting u (j) T ≡ 0 (in an eclipse arc) or using Eqs. (40), (41), and (44) to express the control in terms of state and costate components (along a light arc), until condition (21) is met or t = t f ; if t = t f , then stop the numerical integration, otherwise use the appropriate matching condition (either (47) or (48)) to find the initial value of the adjoint vector for the subsequent arc and repeat steps 5(a) and 5(b); Step 6. Evaluate the violations of the boundary conditions (17) and the necessary conditions (31)- (32). If these violations do not exceed a prescribed threshold, then convergence is declared, otherwise Steps 3 through 6 are repeated.
It is worth remarking that when an eclipse arc is entered, at Step 5(b) all the state and costate components are available in closed form, and numerical integration can be avoided. Details on the costate along coast arcs are provided in Sect. 5. Moreover, the transition time from an eclipse arc to a light arc can be detected in a straightforward and very accurate way. Omitting superscript (j), along coast arcs Eq. (21) can be rewritten as In fact, only x 6,fin and S are time-varying along eclipse arcs. The numerical solution of the preceding equation, in conjunction with the definition of x 6 , the Kepler's equation, and the relation between true anomaly and eccentric anomaly, allows identifying the transition time t j from a light arc to an eclipse arc, without any need of detecting it during the numerical integration process.
The preceding algorithmic steps point out that the optimization process is reduced to identifying the values of some unknown quantities, which form the parameter set, such that all the necessary conditions are fulfilled. This general approach characterizes most indirect optimization algorithms. Using the two relations (47) and (48) and the multipoint conditions (20) and (21), the multiple-arc trajectory optimization problem of interest can be solved through sequential integration of the state and costate equations, while using Eqs. (40), (41), and (44) for the optimal control (in light arcs). Enforcement of Eqs. (27)- (29), (20), and (21) is guaranteed during the integration process. Although in general the parameter set is problem-dependent, for the multiple-arc problem at hand it includes at most the time of flight t f and the initial values of some state components and some (or all the) components of 0 . No intermediate value of any variable belongs to the parameter set. The iterated selection mentioned at Steps 3 and 4 is a delicate operation, which strictly depends on the specific indirect algorithm and regards the parameter set. As an example, the indirect heuristic method [37] employs a heuristic technique (e.g. the particle swarm method) to select the unknown parameters, with the final aim of enforcing all the necessary conditions for optimality, while satisfying the boundary conditions of the problem of interest.
In conclusion, the sequential solution of the multipoint conditions (20) and (21) and the multipoint necessary conditions for optimality (27)-(29) allows identifying a reduced parameter set, which essentially includes the same quantities needed for the solution of a single-arc trajectory optimization problem.

Numerical Example
The methodology and the algorithmic steps described in the previous section are applied to an illustrative numerical example. A low-thrust Earth orbit transfer problem is considered. The initial and final orbits are equatorial and circular, with radii equal to 6778 km and 20,000 km, respectively. Therefore, the boundary condition vector is (cf. Eq. (17)) (49) x 6,fin , t j = 0 where p 0 = 7000 km and p f = 20000 km . Because the two terminal orbits and the transfer path are coplanar, only the state components x 1 , x 2 , x 3 , x 6 , and x 7 (together with the corresponding adjoint variables 1 , 2 , 3 , 6 , and 7 ) are needed in the numerical solution process. This means that the remaining state and costate components are identically zero. The following parameters are assumed for the low-thrust propulsion system: It is worth remarking that the value assumed for u (max) T is beyond the current technological capabilities (although it might become feasible in one or two decades), and was chosen for illustrative purposes.
The boundary condition (30) yields where j j=1,…,7 represent the components of . Therefore, the initial values of the adjoint variables are unknown, with the only exception of 6 . No relation can be identified among the initial values because each of them is expressed in terms of a different component of . Hence, with reference to Step 2 of the numerical solution process described in Sect. 3.3, the minimal set of unknown initial values of the state and costate components is x 6,0 , 1,0 , 2,0 , 3,0 , 7,0 .
In the numerical solution process, canonical units are used: the distance unit (DU) equals 100,000 km, whereas the time unit is such that = 1 DU 3 TU 2 . The indirect heuristic method (IHM) [37], used to solve the problem at hand, assumes the

Minimum-Fuel Trajectories
In most mission scenarios of practical interest, spacecraft are equipped with a finitethrust propulsion system. In these dynamical contexts, the crucial objective consists in minimizing propellant consumption. Previous (and rather extensive) researches proved that minimum-fuel trajectories include relatively short finite-thrust arcs and long-duration coast intervals [30,40]. This section considers the problem of minimizing the propellant consumption for performing an orbit transfer between two specified (initial and final) orbits, while using the modified equinoctial elements to describe the spacecraft dynamics.

Statement of the Problem
The spacecraft of interest is governed by the state equations (11) and is subject to some (problem-dependent) boundary conditions of the form (17). These conditions usually include the relations that define the initial and final orbits. The initial time t 0 is assumed specified and is set to 0. The objective function J to minimize is the propellant mass, and this is equivalent to maximizing the final mass ratio x 7,f . Thus, for the problem at hand the objective is Therefore, the problem consists in finding the optimal control u that minimizes the objective function (54), while holding the state equations (11) and the boundary conditions (17).

Necessary Conditions for Optimality
In order to derive the necessary conditions for optimality, a Hamiltonian function H and a function Φ are defined as (54) J = −x 7,f and the extended objective function J is introduced, Then, the first differential dJ can be obtained by using the analytical steps described in Ref. 22 (similar to those illustrated in Appendix for the more complicated problem of Sect. 3), and is The first differential must vanish at an extremal [22], for arbitrary values of d 0 , d f , dt f , , and , and this implies that the following necessary conditions must hold at an optimal solution: The last condition, which implies stationarity of H with respect to , can be replaced again by the more general Pontryagin minimum principle [32], where the subscript * denotes the optimal value of the corresponding variable. Equa-

3
Journal of Optimization Theory and Applications (2021) 191:545-574 The latter relation, in conjunction with the final condition (66), implies that 7 cannot be positive at all times. Moreover, using the Pontryagin minimum principle and Eq. (69), the optimal value of u T is This means that a minimum-fuel path includes powered phases (where the maximum available thrust is used) and coast arcs. In the previous relation, S is referred to as the switching function, because it determines the switching times between the two types of arcs that compose the optimal trajectory. Equation (71) is deduced under the assumption of neglecting singular arcs, associated with the condition S ≡ 0 over a time interval of finite duration. The existence of similar arcs can be investigated using singular optimal control theory [4]. Unlike the preceding problem, where minimum-time paths including eclipse arcs were sought, minimum-fuel space trajectories do not require formulating a multiple-arc optimization problem. The adjoint vector is continuous for the entire time of flight, and in fact no matching condition for is derived from the necessary conditions for optimality. It is worth stressing that the existence of coast arcs and powered phases along minimum-fuel space trajectories was already proven using different representations for the dynamical state (e.g., Cartesian or spherical coordinates [14,18,20,30,32,34,40]). Therefore, the previous analytical developments represent an alternative derivation leading to an expected result.
In the end, the necessary conditions for optimality (58)-(61) and (63), in conjunction with the state equations (11) and the boundary conditions (17), allow converting the original optimal control problem into a TPBVP, where the unknowns are the state , the control , the final time t f , and the adjoint variables and . An indirect algorithm applied to the problem at hand can proceed using Steps 1-4 and 6 of Sect. 3.2 (with the necessary conditions of the present problem in place of the respective ones found for problem (a)). Instead, Step 5 simplifies in the following fashion: Step 5. Until the current time t ≤ t f , integrate numerically the state and costate equations (11) and (61), while using Eqs.(67), (68), and (71) to express the control in terms of state and costate components.

Closed Form of the Costate Along Optimal Coast Arcs
The preceding two sections address two different trajectory optimization problems that involve multiple coast arcs, where no propulsion is used. Under the assumption of neglecting orbit perturbations, along a coast arc the space vehicle travels a Keplerian trajectory, either an elliptic or a parabolic or a hyperbolic arc. This is apparent also from inspection of Eqs. (2)- (6). In fact, if a r = a = a h = 0 , then . Only x 6 varies along a Keplerian arc (cf. Eq. (7)), due to the true anomaly f (whereas Ω and remain constant). For all the three types of Keplerian paths, f can be found in terms of the current time by solving a transcendental equation (e.g., the Kepler's equation for ellipses and the Barker's equation for parabolas) [12,41]. The time variation of x 6 can be obtained as a result. This section is concerned with the derivation of the closed-form expressions for the adjoint variables j (j = 1, … , 7) along coast arcs. Availability of analytical relations allows precise evaluation of the time histories for all the costate components, and leads to writing the switching function S (cf. Eq. (71)) along coast arcs of minimum-fuel paths as a closed-form expression. As a favorable consequence, the accurate detection of the switching times from coast to thrust arcs, which was recognized as a challenging task in the scientific literature [34], can be facilitated (e.g., using Newton or even higher-order root-finding methods [21], applied to the switching function S).
The Hamiltonian H and the adjoint (or costate) equations (33) and (61) are identical for the two problems addressed in Sects. 3 and 4. The only difference resides in the introduction of multiple arcs for minimum-time problems that include eclipse arcs. However, this has no effect on H and the form of the costate equations. Therefore, for both problems, along coast arcs the Hamiltonian function is Equations (61) and (72) yield the following scalar adjoint equations: As a first result, Eq. (76) proves that the adjoint variables 4 , 5 , and 7 are constant along coast arcs. Therefore, only the remaining components must be integrated. The use of equinoctial elements has thus the remarkable advantage of reducing to 5 the number of state and costate variables that are time-varying along coast arcs, i.e. x 6 , 1 , 2 , 3 , and 6 . This desirable circumstance is not encountered when alternative representations for the state are employed [18,34].
For the purpose of obtaining a closed-form solution to the equation system (73)-(75) and (77), these relations are rewritten with the use of x 6 as the independent variable, in place of the time t. This is possible because x 6 is a strictly monotonic variable for coast arcs of elliptic, parabolic, or hyperbolic type, with positive time derivative ̇x 6 (cf. Eq. (7) with a h = 0 ). Hence, using Eq. (7) (with a h = 0 ), one obtains where the symbol ′ denotes the derivative with respect to x 6 . Equation (81) can be solved by separation of variables. In fact, Eq. (81) can be rewritten as Because x 2 and x 3 are constant along coast arcs, both sides of Eq. (82) are integrable. After a few straightforward steps, the following closed-form solution is obtained: where subscript in denotes the value of the respective variable when the coast arc begins. Insertion of Eq. (83) into Eqs. (78) through (80) yields the following relations: 1 + x 2 cos x 6,in + x 3 sin x 6,in 1 + x 2 cos x 6 + x 3 sin x 6 2 It is worth remarking that Eqs. (83)-(86) hold for Keplerian (coast) arcs of elliptic, parabolic, and hyperbolic type. In the subsections that follow these three types of conics are distinguished, with the intent of obtaining closed-form solutions for 1 , 2 , and 3 .

Elliptic Arcs
In most orbit transfer problems of practical interest, intermediate coast arcs are of elliptic type. This means that x 2 2 + x 2 3 < 1 (cf. the definitions (1)). In this case, letting c 0 ∶= 6,in 1 + x 2 cos x 6,in + x 3 sin x 6,in 2 and using the definition of (cf. Section 2), Eqs. (84)-(86) admit the following closed-form solutions, obtained with the use of the software Mathematica: The symbols c 1 c 2 , and c 3 denote the three integration constants. Their values can be obtained by evaluating Eqs. (87)-(89) at the initial time of the coast arc.

Parabolic Arcs
Parabolic arcs are rather infrequent to encounter in space trajectory optimization. However, for the sake of completeness, this subsection addresses the derivation of the closed-form solution of Eqs. (84)-(86) along parabolic arcs, for which x 2 2 + x 2 3 = 1. As a preliminary step, because x 2 2 + x 2 3 = 1 , the terms 1 + x 2 cos x 6 + x 3 sin x 6 in Eqs. Again, the symbols c 1 c 2 , and c 3 denote the three integration constants. Their values can be obtained by evaluating Eqs. (91)-(93) at the initial time of the coast arc.

Hyperbolic Arcs
Closed-form expressions for the adjoint variables 1 , 2 , and 3 can be obtained also when the space vehicle travels an optimal hyperbolic coast arc, in which x 2 2 + x 2 3 > 1 . In fact, using the identity arctan (i ) = iarctanh (where is a generic variable and i denotes the imaginary unit), Eqs. (87)-(89) can be rewritten for hyperbolic coast arcs as (90) 1 + x 2 cos x 6 + x 3 sin x 6 = 1 + sin x 6 + , with sin = x 2 and cos = x 3 (91) Journal of Optimization Theory and Applications (2021) 191:545-574 you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.