Abstract
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.
Similar content being viewed by others
1 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 formulation of the problem. A similar approach is due to Graham and Rao [19], who employed adaptive Gaussian quadrature orthogonal collocation, without any use of averaging. Their method finds an initial guess by identifying the optimal path while neglecting eclipsing. Then, a sequence of optimal control problems is defined and solved, by adding a single eclipse period at a time. Alternative (indirect) approaches are based on regularization, applied to the transitions from light to shadow (and vice versa). Ferrier and Epenoy [15] used a smoothing function when the spacecraft travels the region between two cylinders that enclose the eclipse cylinder. They proposed an effective methodology for circumventing the theoretical difficulties inherent to the state-dependent control. Most recently, Mazzini [33] utilized a smooth eclipse function, to restore the regularity conditions, for the purpose of applying the Pontryagin minimum principle, whereas Taheri and Junkins [50] proposed the hyperbolic tangent as a smoothing function. Similarly, Aziz et al. [1] smoothed the eclipse transition with a logistic function and computed the optimal paths through differential dynamic programming, a second-order gradient-based method.
The challenges related to the inclusion of multiple coast arcs, associated with eclipse intervals, are thus apparent from the preceding considerations. However, another class of optimal space trajectories exists that may include multiple coast arcs. In fact, minimum-fuel paths using finite thrust are associated with a set of necessary conditions that admit coast arcs and powered phases [40]. This consolidated property was proven in the 60 s [30] and since then a vast amount of literature was dedicated to investigating minimum-fuel space trajectories. A very interesting recent contribution is due to Taheri and Junkins [49], who analyzed the relations between impulsive transfers and finite-thrust paths using optimal control theory. They point out the existence of optimal trajectories with a variety of structures (i.e. different numbers and timing of powered phases and thrust arcs). In this context, the switching function, which depends on the state and the costate variables, plays a major role. In fact, the time evolution of this function determines the sequence of powered phases and coast arcs along the optimal path. For a space vehicle that orbits a single celestial body, orbital motion along thrust intervals must be obtained through numerical integration, whereas its trajectory is Keplerian along coast arcs (under the assumption of neglecting orbit perturbations). This means that the dynamical state is integrable. In the 60 s, Lawden [30] derived the first two-dimensional closed-form also for the costate along optimal coast arcs. Alternative closed-form solutions in the plane of the Keplerian arc are due to Hempel [20], Eckenwiler [14], and Lion and Handlesman [31]. Later, Glandorf [18] derived the general three-dimensional closed-form costate using Cartesian coordinates. More recently, Pan et al. [34] provided the closed-form costate along optimal coast arcs employing spherical coordinates and emphasized its utility for accurate determination of the switching times from coast to thrust.
The analytical study that follows is concerned with both types of trajectory optimization problems that admit multiple coast arcs: (a) minimum-time low-thrust paths with eclipse intervals and (b) minimum-fuel trajectories that employ finite thrust. Modified equinoctial elements are selected as the variables that describe the orbital motion of the space vehicle of interest, subject to the gravitational attraction of a single celestial body. This choice is related to three remarkable properties. 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.
2 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) \(\Omega\), argument of periapse \(\omega\), 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, with the only exception of \(i = \pi\) (condition that is unlikely to encounter, because equatorial retrograde orbits are rather impractical). If \(\vartheta : = 1 + x_{2} \cos x_{6} + x_{3} \sin x_{6}\), the instantaneous radius is \(r = {{x_{1} } \mathord{\left/ {\vphantom {{x_{1} } \vartheta }} \right. \kern-\nulldelimiterspace} \vartheta }\). The classical orbit elements can be retrieved by inverting Eq. (1). The spacecraft position can be written in terms of a, e, i, \(\Omega\), \(\omega\), and f [41] or can be computed directly from the equinoctial elements [19].
2.1 Equations of Motion
The dynamical evolution of the modified equinoctial elements is governed by the respective Gauss equations [3],
where \(\mu\) is the gravitational parameter of the attracting body. The terms \(\left\{ {a_{r} ,a_{\theta } ,a_{h} } \right\}\) are the components of the non-Keplerian acceleration a in the local vertical local horizontal (LVLH) frame aligned with \(\left( {\hat{r},\hat{\theta },\hat{h}} \right)\), where unit vector \(\hat{r}\) is directed toward the instantaneous position vector r (taken from the center of the attracting body), whereas \(\hat{h}\) 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} } \mathord{\left/ {\vphantom {{a_{T} = u_{T} m_{0} } m}} \right. \kern-\nulldelimiterspace} m} = {{u_{T} } \mathord{\left/ {\vphantom {{u_{T} } {x_{7} }}} \right. \kern-\nulldelimiterspace} {x_{7} }}\) and is constrained to the interval \(0 \le a_{T} \le a_{T}^{{\left( {max} \right)}}\), where \(a_{T}^{{\left( {max} \right)}} = {{u_{T}^{{\left( {max} \right)}} } \mathord{\left/ {\vphantom {{u_{T}^{{\left( {max} \right)}} } {x_{7} }}} \right. \kern-\nulldelimiterspace} {x_{7} }}\). The thrust acceleration \(\left( {{{\mathbf{T}} \mathord{\left/ {\vphantom {{\mathbf{T}} m}} \right. \kern-\nulldelimiterspace} m}} \right)\) is assumed to be the only non-Keplerian contribution in Eqs. (2)-(7), thus \({\mathbf{a}} = {{\mathbf{T}} \mathord{\left/ {\vphantom {{\mathbf{T}} m}} \right. \kern-\nulldelimiterspace} m} = {{\mathbf{T}} \mathord{\left/ {\vphantom {{\mathbf{T}} {\left( {m_{0} x_{7} } \right)}}} \right. \kern-\nulldelimiterspace} {\left( {m_{0} x_{7} } \right)}} = {{{\mathbf{u}}_{T} } \mathord{\left/ {\vphantom {{{\mathbf{u}}_{T} } {x_{7} }}} \right. \kern-\nulldelimiterspace} {x_{7} }}\), where \({\mathbf{u}}_{T}\) \(\left( { = {{\mathbf{T}} \mathord{\left/ {\vphantom {{\mathbf{T}} {m_{0} }}} \right. \kern-\nulldelimiterspace} {m_{0} }}} \right)\) has magnitude constrained to the interval \(\left[ {0,u_{T}^{{\left( {max} \right)}} } \right]\). The thrust direction is identified by means of the two thrust angles \(\alpha\) \(\left( { - \pi \le \alpha < \pi } \right)\) and \(\beta\) \(\left( { - {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-\nulldelimiterspace} 2} \le \beta \le {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-\nulldelimiterspace} 2}} \right)\),
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
2.2 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 \(\left( {\hat{c}_{1} ,\hat{c}_{2} ,\hat{c}_{3} } \right)\). Vector \(\hat{c}_{1}\) identifies the vernal axis, which corresponds to the intersection of the Earth equatorial plane with the ecliptic plane, whereas \(\hat{c}_{3}\) points toward the Earth rotation axis. In the ECI-frame, 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 \({\text{c}}_{\eta } : = \cos \eta\) and \({\text{s}}_{\eta } : = \sin \eta\) (\(\eta\) 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, \(\varepsilon \, \left( { = 23.4{\text{ deg}}} \right)\) is the ecliptic obliquity, whereas \(\theta_{S}\) identifies the instantaneous position of the Sun. If \(\theta_{S0}\) denotes the value of \(\theta_{S}\) at the initial time \(t_{0}\) (set to 0) and \(\omega_{S} \, \left( { = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } {T_{S} }}} \right. \kern-\nulldelimiterspace} {T_{S} }}; \, T_{S} = {\text{1 year}}} \right)\) is the (constant) angular rate of the Sun in its motion relative to Earth, then \(\theta_{S} = \omega_{S} t + \theta_{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} \gg R_{E}\), \(\phi_{2} \simeq {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-\nulldelimiterspace} 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.
3 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.
3.1 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 \({{\varvec{\upzeta}}}\) denotes a generic column vector, which depends on the final time \(t_{f}\) and on the initial and final state \({\mathbf{x}}_{0}\) and \({\mathbf{x}}_{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 \({{\varvec{\upzeta}}}\) includes 10 components of the form \(\left\{ {x_{j,0} - \overline{x}_{j,0} } \right\}_{j = 1, \ldots ,5}\) and \(\left\{ {x_{j,f} - \overline{x}_{j,f} } \right\}_{j = 1, \ldots ,5}\), 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_{T}^{\left( j \right)} \equiv 0\)), where inequality (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
The symbol \({\mathbf{x}}_{fin}^{\left( j \right)}\) 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 \(\psi^{\left( j \right)}\) depends explicitly on \(t_{j}\) through the variable \(\theta_{S}\), evaluated at \(t_{j}\) \(\left( {{\text{i}}{\text{.e}}{. }\theta_{S,j} = \omega_{S} t_{j} + \theta_{S0} } \right)\).
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].
3.2 Necessary Conditions for Optimality
In order to derive the necessary conditions for optimality, the extended objective function \(\overline{J}\) is defined,
where \({{\varvec{\upsigma}}}\), \({{\varvec{\upupsilon}}}_{j}\), and \(\xi_{j}\) are time-independent adjoint variables conjugate to the multipoint conditions (17), (20), and (21), whereas \({{\varvec{\uplambda}}}^{\left( j \right)}\) is the time-varying costate vector associated with the state equations (19) and \(t_{N} \equiv t_{f}\). N Hamiltonian functions \(H^{\left( j \right)} \, \left( {j = 1, \ldots ,N} \right)\) and a function \(\Phi\) are introduced,
and the extended objective function \(\overline{J}\) is rewritten as
The first differential \(d\overline{J}\) 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 \(\delta\) 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 \(\left\{ {d{\mathbf{x}}_{ini}^{{\left( {j + 1} \right)}} ,{\mathbf{x}}_{fin}^{\left( j \right)} ,dt_{j} } \right\}_{j = 1, \ldots ,N - 1}\), \(d{\mathbf{x}}_{0}\), \(d{\mathbf{x}}_{f}\), \(dt_{f}\), \(\left\{ {\delta {\mathbf{x}}^{\left( j \right)} ,\delta {\mathbf{u}}^{\left( j \right)} } \right\}_{j = 1, \ldots ,N}\), and this implies that the following necessary conditions for optimality must hold:
The last condition, which implies stationarity of \(H^{\left( j \right)}\) with respect to \({\mathbf{u}}^{\left( j \right)}\), can be replaced by the more general Pontryagin minimum principle [32],
where \(H^{\left( j \right)}\) is defined in Eq. (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 \({{\varvec{\upzeta}}}\) 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 \(\left\{ {t_{1} , \ldots ,t_{N - 1} ,t_{N} \equiv t_{f} } \right\}\). 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.
Using Eqs. (2)–(9) and (23), H can be rewritten as
where \({{\varvec{\upchi}}}_{r}\), \({{\varvec{\upchi}}}_{\theta }\), \({{\varvec{\upchi}}}_{h}\), and \({{\varvec{\upchi}}}\) (whose expressions are not reported for the sake of brevity) depend on \({\mathbf{z}}^{\left( j \right)} : = \left[ {\begin{array}{*{20}c} {x_{1}^{\left( j \right)} } & {x_{2}^{\left( j \right)} } & {x_{3}^{\left( j \right)} } & {x_{4}^{\left( j \right)} } & {x_{5}^{\left( j \right)} } & {x_{6}^{\left( j \right)} } \\ \end{array} } \right]^{T}\). The adjoint variable \(\lambda_{7}^{\left( j \right)}\), which is the seventh component of \({{\varvec{\uplambda}}}^{\left( j \right)}\), deserves special attention. Due to Eqs. (33) and (36), the costate equation for \(\lambda_{7}^{\left( j \right)}\) is
Moreover, because \(\psi^{\left( j \right)}\) is independent of \(x_{7}^{\left( j \right)}\) (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} } \mathord{\left/ {\vphantom {{u_{T} } {x_{7} }}} \right. \kern-\nulldelimiterspace} {x_{7} }} \ge 0\), the thrust angles that minimize \(H^{\left( j \right)}\) are given by
Using these expressions for \(\alpha^{\left( j \right)}\) and \(\beta^{\left( j \right)}\), Eqs. (36) and (37) become
The latter relation, in conjunction with the final condition (39) and the continuity condition (38), implies that \(\lambda_{7}^{\left( j \right)}\) is nonnegative at all times. Moreover, because \(\lambda_{7}^{\left( j \right)} \ge 0\), applying the Pontryagin minimum principle to Eq. (42) allows obtaining the optimal value of \(u_{T}\),
It is worth remarking that Eqs. (40), (41), and (44) are meaningful only for light arcs, where the optimal thrust magnitude corresponds to the maximum value \(u_{T}^{{\left( {max} \right)}} \, \left( { = {{T_{max} } \mathord{\left/ {\vphantom {{T_{max} } {m_{0} }}} \right. \kern-\nulldelimiterspace} {m_{0} }}} \right)\), as prescribed by Eq. (44).
In the end, the necessary conditions for optimality (27)-(33) and (35), in conjunction with the state equations (19) and the multipoint conditions (20)-(21), allow converting the original optimal control problem into a two-point boundary value problem (TPBVP), where the unknowns are the state \({\mathbf{x}}^{\left( j \right)} \, \left( {j = 1, \ldots ,N} \right)\), the control \({\mathbf{u}}^{\left( j \right)} \, \left( {j = 1, \ldots ,N} \right)\), the intermediate and final times \(\left\{ {t_{1} , \ldots ,t_{N - 1} ,t_{N} \equiv t_{f} } \right\}\), and the adjoint variables \({{\varvec{\uplambda}}}^{\left( j \right)} \, \left( {j = 1, \ldots ,N} \right)\), \({{\varvec{\upsigma}}}\), \({{\varvec{\upupsilon}}}_{k}\), and \(\xi_{k}\) \(\left( {k = 1, \ldots ,N - 1} \right)\).
3.3 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
Moreover, Eq. (29) can be solved for \(\xi_{j}\), using also Eq. (42) to express the Hamiltonian functions \(H_{ini}^{{\left( {j + 1} \right)}} {\text{ and }}H_{fin}^{\left( j \right)}\),
where \({{\varvec{\upchi}}}_{r}^{{}} : = {{\varvec{\upchi}}}_{r,ini}^{{\left( {j + 1} \right)}} \equiv {{\varvec{\upchi}}}_{r,fin}^{\left( j \right)}\), \({{\varvec{\upchi}}}_{\theta }^{{}} : = {{\varvec{\upchi}}}_{\theta ,ini}^{{\left( {j + 1} \right)}} \equiv {{\varvec{\upchi}}}_{\theta ,fin}^{\left( j \right)}\), \({{\varvec{\upchi}}}_{h}^{{}} : = {{\varvec{\upchi}}}_{h,ini}^{{\left( {j + 1} \right)}} \equiv {{\varvec{\upchi}}}_{h,fin}^{\left( j \right)}\), and \({{\varvec{\upchi}}}: = {{\varvec{\upchi}}}_{ini}^{{\left( {j + 1} \right)}} \equiv {{\varvec{\upchi}}}_{fin}^{\left( j \right)}\), because the state is continuous across adjacent arcs. Two cases can occur, with regard to Eq. (46): (a) transition from a light arc to an eclipse arc, which means that \(u_{T,j}^{{\left( {max} \right)}} = u_{T}^{{\left( {max} \right)}}\) and \(u_{T,j + 1}^{{\left( {max} \right)}} = 0\), and (b) transition from an eclipse arc to a light arc, implying \(u_{T,j}^{{\left( {max} \right)}} = 0\) and \(u_{T,j + 1}^{{\left( {max} \right)}} = u_{T}^{{\left( {max} \right)}}\).
Case (a). Transition from a light arc to an eclipse arc. Insertion of Eq. (46) (with \(u_{T,j}^{{\left( {max} \right)}} = u_{T}^{{\left( {max} \right)}}\) and \(u_{T,j + 1}^{{\left( {max} \right)}} = 0\)) into Eq. (45) leads to obtaining the following relation:
where I denotes the identity matrix, with dimension appropriate to the context. Equation (47) can be regarded as a linear system with \({{\varvec{\uplambda}}}_{ini}^{{\left( {j + 1} \right)}}\) as the unknown. Therefore, \({{\varvec{\uplambda}}}_{ini}^{{\left( {j + 1} \right)}}\) can be obtained using Eq. (47), if all the other quantities (including \({{\varvec{\uplambda}}}_{fin}^{\left( j \right)}\)) are known.
Case (b). Transition from an eclipse arc to a light arc. Insertion of Eq. (46) (with \(u_{T,j}^{{\left( {max} \right)}} = 0\) and \(u_{T,j + 1}^{{\left( {max} \right)}} = u_{T}^{{\left( {max} \right)}}\)) into Eq. (45) leads to obtaining the following vector equation:
Squaring both sides of Eq. (48) leads to obtaining a system of quadratic equations with the components of \({{\varvec{\uplambda}}}_{ini}^{{\left( {j + 1} \right)}}\) as the unknowns. The entire set of solutions of this system can be obtained numerically (e.g., in MATLAB using the embedded function vpasolve). Then, among the real solutions, only those that preserve the sign consistency of Eq. (48) are acceptable. If more than a single acceptable solution for \({{\varvec{\uplambda}}}_{ini}^{{\left( {j + 1} \right)}}\) exists, then the one that minimizes \(\left| {{{\varvec{\uplambda}}}_{ini}^{{\left( {j + 1} \right)}} - {{\varvec{\uplambda}}}_{fin}^{\left( j \right)} } \right|\) is selected. This choice is consistent with a straightforward property that regards \(\left\{ {{{\varvec{\uplambda}}}_{ini}^{{\left( {j + 1} \right)}} ,{{\varvec{\uplambda}}}_{fin}^{\left( j \right)} } \right\}\). In fact, as the eclipse arc duration reduces or as \(u_{T}^{{\left( {max} \right)}}\) tends to zero, \({{\varvec{\uplambda}}}_{ini}^{{\left( {j + 1} \right)}}\) must tend to \({{\varvec{\uplambda}}}_{fin}^{\left( j \right)}\), 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 \({{\varvec{\uplambda}}}\) at times \(t_{j} \, \left( {j = 1, \ldots ,N - 1} \right)\). 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 \({{\varvec{\uplambda}}}\)) by eliminating some components of \({{\varvec{\upsigma}}}\) from Eq. (30); this leads to identifying the minimal set of unknown initial values for the components of \({\mathbf{x}}_{0}\) and \({{\varvec{\uplambda}}}_{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 \le 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_{T}^{\left( j \right)} \equiv 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 \(\theta_{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 \({{\varvec{\uplambda}}}_{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.
3.4 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))
where \(p_{0} = 7000{\text{ km}}\) and \(p_{f} = 20000{\text{ 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 \(\lambda_{1}\), \(\lambda_{2}\), \(\lambda_{3}\), \(\lambda_{6}\), and \(\lambda_{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_{T}^{{\left( {max} \right)}}\) 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 \(\left\{ {\sigma_{j} } \right\}_{j = 1, \ldots ,7}\) represent the components of \({{\varvec{\upsigma}}}\). Therefore, the initial values of the adjoint variables are unknown, with the only exception of \(\lambda_{6}\). No relation can be identified among the initial values because each of them is expressed in terms of a different component of \({{\varvec{\upsigma}}}\). 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 \(\left\{ {x_{6,0} ,\lambda_{1,0} ,\lambda_{2,0} ,\lambda_{3,0} ,\lambda_{7,0} } \right\}\).
In the numerical solution process, canonical units are used: the distance unit (DU) equals 100,000 km, whereas the time unit is such that \(\mu = 1 \, {{{\text{DU}}^{{3}} } \mathord{\left/ {\vphantom {{{\text{DU}}^{{3}} } {{\text{TU}}^{{2}} }}} \right. \kern-\nulldelimiterspace} {{\text{TU}}^{{2}} }}\). The indirect heuristic method (IHM) [37], used to solve the problem at hand, assumes the preceding unknown initial values and the time of flight as the parameter set. A thorough description of IHM can be found in Ref. 37. Further details on the numerical solution process are omitted for the sake of brevity. The value of \(\theta_{S0}\) (cf. Sect. 2.2) is set to 10 deg. The minimum time of flight equals 5.1 days, whereas the boundary conditions do not exceed \(10^{ - 3}\). The magnitude of \({{\varvec{\upzeta}}}\) (cf. Eq. (50)) is \(7.4 \cdot 10^{ - 4}\). Figures 1 and 2 depict the time histories of the semilatus rectum (component \(x_{1}\)) and the orbit eccentricity (retrieved from \(x_{2}\) and \(x_{3}\)). The horizontal segments correspond to the eclipse arcs, where no thrust is employed and the orbit elements remain constant as a result. Figures 3 and 4 portray the time histories of the adjoint variables \(\lambda_{2}\) and \(\lambda_{3}\). Inspection of these figures reveals that at the transition points \(\lambda_{2}\) and \(\lambda_{3}\) are subject to discontinuities, shown with greater detail in the insets. Similar time behaviors characterize the time histories of \(\lambda_{1}\) and \(\lambda_{6}\) (while \(\lambda_{7}\) is continuous, cf. Eq. (38)). The numerical results found for this illustrative example definitely corroborate the effectiveness of the solution methodology described in Sect. 3.3.
Time history of the semilatus rectum (with zoom in the inset)
Time history of the eccentricity (with zoom in the inset)
Time history of the adjoint variable \(\lambda_{2}\) (with zoom in the inset)
Time history of the adjoint variable \(\lambda_{3}\) (with zoom in the inset)
4 Minimum-Fuel Trajectories
In most mission scenarios of practical interest, spacecraft are equipped with a finite-thrust 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.
4.1 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).
4.2 Necessary Conditions for Optimality
In order to derive the necessary conditions for optimality, a Hamiltonian function \(H\) and a function \(\Phi\) are defined as
and the extended objective function \(\overline{J}\) is introduced,
Then, the first differential \(d\overline{J}\) 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{\mathbf{x}}_{0}\), \(d{\mathbf{x}}_{f}\), \(dt_{f}\), \(\delta {\mathbf{x}},{\text{ and }}\delta {\mathbf{u}}\), 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 \({\mathbf{u}}\), can be replaced again by the more general Pontryagin minimum principle [32],
where the subscript * denotes the optimal value of the corresponding variable. Equation (61) is the adjoint equation for the costate vector, accompanied by the related boundary conditions (58) and (59). Because \({{\varvec{\upzeta}}}\) appears in Eqs. (58) and (59), their explicit form is problem-dependent, unlike what occurs for the adjoint equations (61). Furthermore, Eq. (60) holds for the final value of the Hamiltonian function.
Using Eqs. (2)–(9) and (63), H can be rewritten as
where z collects components \(x_{1}\) through \(x_{6}\) of the state, i.e. \({\mathbf{z}}: = \left[ {\begin{array}{*{20}c} {x_{1} } & {x_{2} } & {x_{3} } & {x_{4} } & {x_{5} } & {x_{6} } \\ \end{array} } \right]^{T}\); the terms \(H_{r}\), \(H_{\theta }\), \(H_{h}\), and \(H_{0}\), whose expressions are not reported for the sake of brevity, depend on both z and \({{\varvec{\uplambda}}}\). Due to Eq. (64), the adjoint equation for \(\lambda_{7}\) is
Moreover, because the final mass is unspecified (and in fact is to be minimized), the boundary conditions (17) are independent of \(x_{7,f}\). As a result, Eq. (59) yields
The minimum principle (63) allows expressing the optimal control in terms of the state and costate variables. With reference to Eq. (64), the first three terms in square parentheses can be regarded as a dot product. Thus, since \({{u_{T} } \mathord{\left/ {\vphantom {{u_{T} } {x_{7} }}} \right. \kern-\nulldelimiterspace} {x_{7} }} \ge 0\), the thrust angles that minimize H are given by
Using these expressions for \(\alpha\) and \(\beta\), Eqs. (64) and (65) become
The latter relation, in conjunction with the final condition (66), implies that \(\lambda_{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 \equiv 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 \({{\varvec{\uplambda}}}\) is continuous for the entire time of flight, and in fact no matching condition for \({{\varvec{\uplambda}}}\) 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 \({\mathbf{x}}\), the control \({\mathbf{u}}\), the final time \(t_{f}\), and the adjoint variables \({{\varvec{\uplambda}}}\) and \({{\varvec{\upsigma}}}\).
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 \le 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.
5 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_{\theta } = a_{h} = 0\), then \(\dot{x}_{1} = \dot{x}_{2} = \dot{x}_{3} = \dot{x}_{4} = \dot{x}_{5} = 0\). Only \(x_{6}\) varies along a Keplerian arc (cf. Eq. (7)), due to the true anomaly f (whereas \(\Omega\) and \(\omega\) 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 \(\lambda_{j} \, \left( {j = 1, \ldots ,7} \right)\) 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 \(\lambda_{4}\), \(\lambda_{5}\), and \(\lambda_{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} , \, \lambda_{1} , \, \lambda_{2} , \, \lambda_{3} ,{\text{ and }}\lambda_{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 \(\dot{x}_{6}\) (cf. Eq. (7) with \(a_{h} = 0\)). Hence, using Eq. (7) (with \(a_{h} = 0\)), one obtains
where the symbol \(^{\prime }\) 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:
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 \(\lambda_{1}\), \(\lambda_{2}\), and \(\lambda_{3}\).
5.1 Elliptic Arcs
In most orbit transfer problems of practical interest, intermediate coast arcs are of elliptic type. This means that \(x_{2}^{2} + x_{3}^{2} < 1\) (cf. the definitions (1)). In this case, letting \(c_{0} : = \lambda_{6,in} \left( {1 + x_{2} \cos x_{6,in} + x_{3} \sin x_{6,in} } \right)^{2}\) and using the definition of \(\vartheta\) (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.
5.2 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_{3}^{2} = 1\).
As a preliminary step, because \(x_{2}^{2} + x_{3}^{2} = 1\), the terms \(\left( {1 + x_{2} \cos x_{6} + x_{3} \sin x_{6} } \right)\) in Eqs. (84)–(86) are rewritten as
The preceding definitions of \(c_{0}\) and \(\vartheta\) are used again. After insertion of Eq. (90) into Eqs. (84)–(86), the following closed-form solutions are obtained with the use of the software Mathematica:
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.
5.3 Hyperbolic Arcs
Closed-form expressions for the adjoint variables \(\lambda_{1}\), \(\lambda_{2}\), and \(\lambda_{3}\) can be obtained also when the space vehicle travels an optimal hyperbolic coast arc, in which \(x_{2}^{2} + x_{3}^{2} > 1\). In fact, using the identity \(\arctan \left( {i\eta } \right) = i{\text{arctanh}} \eta\) (where \(\eta\) is a generic variable and i denotes the imaginary unit), Eqs. (87)–(89) can be rewritten for hyperbolic coast arcs as
Once more, the symbols \(c_{1}\) \(c_{2}\), and \(c_{3}\) denote the three integration constants. Their values can be obtained by evaluating Eqs. (94)–(96) at the initial time of the coast arc.
6 Concluding Remarks
This paper presents the analytical study of two types of optimal space trajectories that include multiple coast arcs: (a) minimum-time low-thrust paths with eclipse intervals and (b) minimum-fuel trajectories that employ finite thrust. Modified equinoctial elements are used to describe the orbit dynamics.
Problem (a) is formulated as a multiple-arc optimization problem, and all the related necessary conditions for optimality are derived. These form an extended set of conditions, which includes the Pontryagin minimum principle (in the light arcs), the adjoint equations for the costate variables, the related boundary conditions, the transversality relation on the final value of the Hamiltonian, and the multipoint necessary conditions, at the junction times between two arcs. The latter relations are combined, to yield the matching (jump) conditions for the costate variables between two consecutive arcs. This research demonstrates that all the multipoint conditions can be solved sequentially in the numerical solution process. As a result, the parameter set for an indirect algorithm retains the size of the typical set associated with a single-arc optimization problem. No regularization or averaging is required to make tractable and solve the problem. In fact, based on this new multiple-arc formulation and the related analysis, indirect algorithms can easily incorporate the discontinuities of the costate variables, located at the times when the spacecraft transitions from light to shadow (and vice versa). Moreover, the use of equinoctial elements can mitigate the hypersensitivity of the solution on the initial values of the adjoints, which is a common issue when indirect approaches are used. As a result, an indirect method is capable of yielding very accurate numerical results, while avoiding several difficulties of a theoretical or numerical nature, inherent to alternative formulations. An illustrative numerical example demonstrates the effectiveness of the approach proposed in this work.
This research also revisits problem (b) with the use of equinoctial elements, and derives the necessary conditions associated with minimum-fuel trajectories. Unlike minimum-time paths with eclipse intervals, problem (b) is formulated as a single-arc optimization problem. The switching function is introduced, and plays the role of determining the optimal sequence of powered phases and thrust arcs, whose number and timing is unknown a priori. The well-consolidated properties of minimum-fuel trajectories are retrieved in terms of equinoctial elements and related adjoint variables, and the substantial analytical differences between problems (a) and (b) are emphasized.
As a further contribution, this study focuses on the costate along optimal coast arcs. Closed-form expressions, written in terms of elementary functions, are derived for the adjoint variables associated with modified equinoctial elements, along the three major types of Keplerian trajectories. Overall, 9 out of 14 components of the state and the costate turn out to be constant along optimal coast arcs, whereas for the remaining 5 variables closed-form expressions exist. This allows fast and accurate evaluation of the state and costate time histories along coast arcs, and leads to writing the switching function for problem (b) in closed form. This finding can be beneficial for the accurate determination of the switching times from coast to thrust arcs.
These circumstances definitely represent further, unequivocal advantages of using modified equinoctial elements in spacecraft trajectory optimization.
References
Aziz, J., Scheeres, D., Parker, J., Englander, J.: A smoothed eclipse model for solar electric propulsion trajectory optimization. Trans. Jpn. Soc. Aeronaut. Space Sci. 17(2), 181–188 (2019)
Bastante, J.C., Letizia, F., De Bruijn, F., Lubke-Ossenbeck, B.: Application of the sequential gradient restoration algorithm to the solution of the low thrust transfer problem. Aerotecnica Missili Spazio 96(4), 228–236 (2017)
Battin, R.H. An Introduction to the Mathematics and Methods of Astrodynamics, AIAA Education Series, New York, NY, 1987, pp. 448–450, 490–494
Bell, D.J., Jacobson, D.H.: Singular Optimal Control Problems, pp. 1–100. Academic Press, New York, NY (1975)
Betts, J.T.: Very low-thrust trajectory optimization using a direct SQP method. J. Comput. Appl. Math. 120, 27–40 (2000)
Betts, J.T., Erb, S.O.: Optimal low thrust trajectories to the moon. SIAM J. Appl. Dyn. Syst. 2(2), 144–170 (2003)
Betts, J.T.: Optimal low thrust orbit transfers with eclipsing. Optim. Control Appl. Methods 36, 218–240 (2015)
Betts, J.T.: Survey of numerical methods for trajectory optimization. J. Guid. Control. Dyn. 21(2), 193–207 (1988)
Broucke, R.A., Cefola, P.J.: On the equinoctial orbit elements. Celest. Mech. 5, 303–310 (1972)
Cerf, M.. Fast solution of minimum-time low-thrust transfer with eclipses. In: Proceedings of the institution of mechanical engineers part G: Journal of Aerospace Engineering, Vol. 233, No. 7, 2019, pp. 2699–2714
Conway, B.A.: A survey of methods available for the numerical optimization of continuous dynamical systems. J. Optim. Theory Appl. 152, 271–306 (2012)
Curtis H.D.. Orbital Mechanics for Engineering Students, Elsevier, Kidlington, Oxford, U.K., 2014, pp. 147–169, 701–702
Edelbaum, T.N.: The use of high- and low-thrust propulsion in combination for space missions. J. Astronaut. Sci. 9, 58–59 (1962)
Eckenwiler, M.W.: Closed-form Lagrangian multipliers for coast periods of optimum trajectories. AIAA J. 3(6), 1149–1151 (1965)
Ferrier, C., Epenoy, R.: Optimal control for engines with eletro-ionic propulsion under constraint of eclipse. Acta Astronaut. 48(4), 181–192 (2001)
Gao, Y.: Near-optimal very low-thrust earth orbit transfers and guidance schemes. J. Guid. Control. Dyn. 30(2), 529–539 (2007)
Geoffrey, S., Epenoy, R.: Optimal low-thrust transfers with constraints – Generalization of averaging techniques. Acta Astronaut. 48(4), 181–192 (2001)
Glandorf, D.R.: Lagrange multipliers and the state transition matrix for coasting arcs. AIAA J. 7(2), 363–365 (1968)
Graham, K.F., Rao, A.V.: Minimum-time trajectory optimization of low-thrust earth-orbit transfers with eclipsing. J. Guid. Control. Dyn. 53(2), 289–303 (2016)
Hempel, P.R.: Representation of the Lagrangian multipliers for coast periods of optimum trajectories. AIAA J. 4(4), 729–730 (1966)
Hildebrand, F.B.: Introduction to Numerical Analysis, pp. 578–583. Dover, Mineola, NY (1987)
D. G. Hull. Optimal Control Theory for Applications, Springer International Edition, New York, NY, 2003, pp. 95–100, 142–143, 221–223
Kéchichian, J.A.: Optimal low-thrust rendezvous using equinoctial orbit elements. Acta Astronaut. 38(1), 1–14 (1996)
Kéchichian, J.A.: Optimal low-thrust transfer in general circular orbit using analytic averaging of the system dynamics. J. Astronaut. Sci. 57, 369–392 (2009)
Kéchichian, J.A.: Mathematics of attitude-constrained optimal low-thrust orbit transfer. Aerotecnica Missili Spazio 96(4), 180–194 (2017)
Kéchichian, J.A.: Orbit raising with low-thrust tangential acceleration in presence of earth shadow. J. Spacecr. Rocket. 35(4), 516–525 (1998)
Kéchichian, J.A.: Low-thrust inclination control in presence of earth shadow. J. Spacecr. Rocket. 35(4), 526–532 (1998)
Kluever, C.A., Oleson, S.R.: Direct approach for computing near-optimal low-thrust earth-orbit tranfers. J. Spacecr. Rocket. 35(4), 509–515 (1998)
Kluever, C.A.: Using Edelbaum’s method to compute low-thrust transfers with earth-shadow eclipses. J. Guid. Control. Dyn. 34(1), 300–303 (2011)
Lawden, D.F.: Optimal Trajectories for Space Navigation. Butterworths, London, U. K. (1963)
Lion, P.M., Handelsman, M.: Primer vector on fixed-time impulsive trajectories. AIAA J. 6(1), 127–132 (1968)
Longuski, J.M., Guzman, J.J., Prussing, J.E.: Optimal Control with Aerospace Applications. Springer, New York (2014)
Mazzini, L.: Finite thrust orbital transfers. Acta Astronaut. 100, 107–128 (2014)
Pan, B., Lu, P., Chen, Z.: Three-dimensional closed-form costate solutions in optimal coast. Acta Astronaut. 77, 156–166 (2012)
Petropoulos, A.E., Simple control laws for low-thrust orbit transfers. AIAA/AAS Astrodynamics Specialist Conference, Big Sky, MN, 2003
Petropoulos A.E., Low-thrust orbit transfers using candidate Lyapunov functions with a mechanism for coasting. AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Providence, RI, 2004
Pontani, M., Conway, B.A.: Optimal low-thrust orbital maneuvers via indirect swarming method. J. Optim. Theory Appl. 162, 272–292 (2014)
Pontani, M.: Optimal low-thrust trajectories using nonsingular equinoctial orbit elements. Adv. Astronaut. Sci. 170, 443–462 (2020)
Pontani M., Metodi di Ottimizzazione Globale delle Traiettorie Aerospaziali, Ph.D. dissertation, Sapienza Università di Roma, 2008, pp. 10–17
Prussing, J.E.: "Primer Vector Theory and Applications. In: Conway, B. (ed.) Spacecraft Trajectory Optimization, pp. 16–36. Cambridge University Press, New York (2001)
Prussing, J.E., Conway, B.A.: Orbital Mechanics. Oxford University Press, New York (2013)
Rao A.V. A survey of numerical methods for optimal control. Advances in the Astronautical Sciences, Vol. 135, 2010; paper AAS 09–334
Rathsman, P., Kugelberg, J., Bodin, P., Racca, G.D., Foing, B., Stagnaro, L.: SMART-1: Development and lessons learnt. Acta Astronaut. 57(2–8), 455–468 (2005)
Rayman, M.D., Chadbourne, P.A., Culwell, J.S., Williams, S.N.: Mission design for deep space 1: a low-thrust technology validation mission. Acta Astronaut. 45(4–9), 381–388 (1999)
Robbins, H.M.: An Analytical Study of the Impulsive Approximation. AIAA J. 4(8), 1417–1423 (1966)
Ross, I.M., Gong, Q., Sekhavat, P.: Low-thrust, high-accuracy trajectory optimization. J. Guid. Control. Dyn. 30(4), 921–933 (2007)
Singh S.K., Taheri E., Woollands R., Junkins J., Mission design for close-range lunar mapping by quasi-frozen orbits. Proceedings of the 70th International Astronautical Congress, Washington, D.C., 2019; paper IAC-19-C1.1.11
Sutton, G.P., Biblarz, O.: Rocket Propulsion Elements, pp. 660–710. John Wiley and Sons, New York, NY (2001)
Taheri, E., Junkins, J.L.: How many impulses redux. J. Astronaut. Sci. 67, 257–334 (2020)
Taheri, E., Junkins, J.L.: Generic smoothing for optimal bang-off-bang spacecraft maneuvers. J. Guid. Control. Dyn. 41(11), 2467–2472 (2018)
Funding
Open access funding provided by Università degli Studi di Roma La Sapienza within the CRUI-CARE Agreement.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Communicated by: Francesco Topputo.
Appendices
Appendix
First Differential of the Extended Objective Function for Multiple-Arc Problems
This Appendix is concerned with the derivation of the first differential of the extended objective function (25), in the context of the multiple-arc optimization problem addressed in Sect. 3. The dynamical system of interest is governed by state equations of the form (19) and is subject to the boundary conditions (17).
As a first step, the first differential of the scalar term \(d\Phi\) is obtained, using the chain rule,
However, all the terms in the second row of Eq. (97) vanish. In fact, due to Eq. (24),
As a second step, the differential of the remaining terms of (25) is derived, using the general relation for the differential of an integral [22],
where the symbol \(\delta\) denotes the variation, i.e. the time-fixed differential [22]. However, the last term in the second row of Eq. (99) vanishes, because, due to Eq. (19)
In Eq. (99) the following term is integrated by parts:
The relation \(d\eta = \delta \eta + \dot{\eta }dt\) (where \(\eta\) represents a generic variable) [22] is used in Eq. (101), then the resulting expression is inserted into Eq. (99), to yield
Finally, summation of the two right-hand sides of Eqs. (97) and (102) (while holding Eqs. (98) and (100)) leads to obtaining the final expression of the first differential \(d\overline{J}\), i.e. Eq. (26).
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as 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://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Pontani, M. Optimal Space Trajectories with Multiple Coast Arcs Using Modified Equinoctial Elements. J Optim Theory Appl 191, 545–574 (2021). https://doi.org/10.1007/s10957-021-01867-2
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10957-021-01867-2