Minimum-Energy Transfer Optimization between Near-Circular Orbits Using an Approximate Closed-Form Solution

Abstract

This paper explores the use of the averaging method in the optimal control problem related to the multirevolution orbital transfer of a spacecraft with low-thrust capabilities. The regularized equations of motion are expressed using modified equinoctial elements with the eccentric longitude as a fast variable. The control function is represented as a Fourier series relative to the eccentric longitude. The classical averaging technique’s usage results in the averaged trajectory depending only on a limited number of optimization parameters. Moreover, when transferring between near-circular orbits, the averaged motion can be estimated using analytical formulas. As such, the optimal multiorbit flight problem is simplified to nonlinear programming with fewer parameters, thereby accelerating the optimal solution’s derivation. Two practical examples illustrate the technique’s application: orbital transfer near the geostationary orbit and circular orbit raising maneuver. The solutions derived are compared with Pontryagin extremals.

1. Introduction

Many space missions prefer to utilize low-thrust propulsion systems because they have higher specific impulse compared with high-thrust propulsion systems that significantly decreases the quantity of the propellant used. Two commonly used approaches to low-thrust optimal control problems are indirect and direct methods [1,2]. Indirect methods utilize the Pontryagin maximum principle, which transforms the optimal control problem into a two-point boundary value problem [3,4,5]. Using indirect optimization, the global search for fuel-optimal low-thrust transfers from a low Earth orbit to the vicinity of the Moon was performed in the Earth–Moon planar circular restricted three-body problem [6]. Extensive research and classification of low-thrust low-energy transfers between low Earth and lunar orbits using indirect optimal control are thoroughly examined in [7]. On the other hand, direct methods rely on discretization/parametrization and transform the optimal control problem into a nonlinear programming problem [8,9,10] to solve transfer [11], rephasing [12], and rendezvous [13,14,15,16] problems. Both approaches struggle with significant sensitivity to the initial guess, lead to numerous local extrema, and incur high computational costs during the optimization process. To overcome these difficulties, various methods are used [17], including modified continuation techniques [18], search space landscape analysis [19], convexification [20,21], and neural control [22,23].

Both direct and indirect methods often use averaging as a key preliminary step in solving low-thrust optimization problems. The averaged equations provide an overview of the orbit’s mean evolution, excluding short-term oscillations resulting from the thrust or disturbances. Additionally, they are numerically more stable than the system without averaging. The application of averaging allows for a fast and accurate search for optimal orbital transfers [24,25,26]. In [27,28], analytical solutions of minimum fuel and minimum energy transfer problems were achieved for close elliptic coplanar orbits using the Pontryagin principle with the averaging of the Hamiltonian.

Averaging is also applied within the direct optimization framework. According to a study [29], when the thrust (or perturbation) is a function with a constant spectrum, the averaged right-hand side of the Gauss variational equations depends solely on 14 thrust-perturbation Fourier coefficients. This fact was used to define low-thrust control using a finite number of variables and solve an optimal control problem using direct methods [30]. However, the classical orbital elements used in [29] are singular in the instances of circular and equatorial orbits. Additionally, the averaged variational equations in classical orbital elements are not integrable.

The authors present an advancement of the above approach in the current study. By expressing the spacecraft orbital motion in terms of modified equinoctial elements and taking the eccentric longitude as a fast variable, we found that the averaged equations of motion contain only 13 Fourier coefficients. Moreover, a closed-form solution for cases with low eccentricity is found. Now we are interested in the application of this averaging technique with analytical formulas to multiorbit low-thrust transfer optimization problems. The essence of this work is the reduction of the optimal control problem to a nonlinear programming problem with a small number of optimization parameters and its solution.

This paper is organized as follows. Section 2 presents the general optimal control problem that is under investigation. Section 3.1 explains the developed averaging approach. Then, in Section 3.2, the case of transfer between near-circular orbits is formalized, and the corresponding closed-form solution of the averaged equations is presented. Section 4 outlines the process of transforming the optimal control problem into nonlinear programming using the averaging method. In Section 5, two test energy-optimal control problems are defined and solved: orbital transfer near the geostationary orbit (Section 5.1) and circular orbit raising maneuver (Section 5.2). Section 5.3 compares the optimal solutions with the Pontryagin extremals, which were obtained using the Pontryagin maximum principle. The final section draws conclusions from the study and suggests directions for further research on the topic.

2. The General Optimal Control Problem under Study

The perturbed orbital motion of the spacecraft in the vicinity of the center of attraction is considered. The motion of the spacecraft is described by a system of differential equations in terms of the modified equinoctial elements as follows [31]:

$$\begin{array}{cc}\hfill \frac{dp}{dt}& =\frac{2}{\sigma \left(F\right)}\sqrt{\frac{p^3}{\mu }}·f_c,\hfill \\ \hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{d{e}_x}{dt}& =\sqrt{\frac{p}{\mu }}\left\{sinL\left(F\right)·f_r+\left(cosL\left(F\right)+\frac{e_x+cosL\left(F\right)}{\sigma \left(F\right)}\right)·f_c-e_y\frac{i_{x}sinL\left(F\right)-i_{y}cosL\left(F\right)}{\sigma \left(F\right)}·f_n\right\},\hfill \\ \hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{d{e}_y}{dt}& =\sqrt{\frac{p}{\mu }}\left\{-cosL\left(F\right)·f_r+\left(sinL\left(F\right)+\frac{e_y+sinL\left(F\right)}{\sigma \left(F\right)}\right)·f_c+e_x\frac{i_{x}sinL\left(F\right)-i_{y}cosL\left(F\right)}{\sigma \left(F\right)}·f_n\right\},\hfill \\ \hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{d{i}_x}{dt}& =\sqrt{\frac{p}{\mu }}\frac{1+i_x^2+i_y^2}{2\sigma \left(F\right)}cosL\left(F\right)·f_n,\hfill \\ \hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \frac{d{i}_y}{dt}& =\sqrt{\frac{p}{\mu }}\frac{1+i_x^2+i_y^2}{2\sigma \left(F\right)}sinL\left(F\right)·f_n,\hfill \\ \hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \frac{dF}{dt}& ={\varphi }^3\sqrt{\frac{\mu }{p^3}}+\frac{d{e}_x}{dt}cosF-\frac{d{e}_y}{dt}sinF+\hfill \\ & +\sqrt{\frac{p}{\mu }}\left\{-\frac{2}{\varphi }·f_r-b\frac{1+\sigma \left(F\right)}{{\varphi }^2}(e_{y}cosL\left(F\right)-e_{x}sinL\left(F\right))·f_c+\frac{i_{x}sinL\left(F\right)-i_{y}cosL\left(F\right)}{{\varphi }^2}·f_n\right\},\hfill \end{array}$$

(6)

where

$$\begin{array}{cc}\hfill \sigma \left(F\right)& =1+e_{x}cosL\left(F\right)+e_{y}sinL\left(F\right),\hfill \\ \hfill sinL\left(F\right)& =\frac{(1-e_x^{2}b)sinF+e_x{e}_{y}bcosF-e_y}{1-e_{x}cosF-e_{y}sinF},\hfill \\ \hfill cosL\left(F\right)& =\frac{(1-e_y^{2}b)cosF+e_x{e}_{y}bsinF-e_x}{1-e_{x}cosF-e_{y}sinF},\hfill \\ \hfill b& =\frac{1}{1+\varphi },\hfill \\ \hfill \varphi & =\sqrt{1-e_x^2-e_y^2}.\hfill \end{array}$$

Here, p, $e_x$, $e_y$, $i_x$, and $i_y$ are the modified equinoctial elements that are related to the semimajor axis a, eccentricity e, inclination i, right ascension of the ascending node $\mathrm{\Omega }$, and argument of periapsis $\omega$ by the expressions [31]

$$\begin{array}{cc}\hfill p& =a(1-e^2),\hfill \\ \hfill e_x& =ecos(\mathrm{\Omega }+\omega ),\hfill \\ \hfill e_y& =esin(\mathrm{\Omega }+\omega ),\hfill \\ \hfill i_x& =tan(i/2)cos\mathrm{\Omega },\hfill \\ \hfill i_y& =tan(i/2)sin\mathrm{\Omega }.\hfill \end{array}$$

(7)

The variable F is the eccentric longitude; it relates to $\mathrm{\Omega }$, $\omega$, and the eccentric anomaly E by the expression $F=\mathrm{\Omega }+\omega +E$. The perturbation variables $f_r$, $f_c$, and $f_n$ are the projections of the perturbation acceleration $\mathbf{f}$ on the axes of the local-vertical/local-horizontal frame; $f_r$ is the radial component; $f_c$ is the component orthogonal to the radius vector of the spacecraft; and $f_n$ is the component normal to the orbital plane. The parameter $\mu$ is the gravitational parameter of the attracting body.

The equations above can be written in a compact form:

$$\begin{array}{c}\hfill \dot{\mathbf{x}}=\mathbf{g}(\mathbf{x},F,\mathbf{f}),\\ \hfill \end{array}$$

(8)

$$\begin{array}{c}\hfill \dot{F}=h(\mathbf{x},F,\mathbf{f}),\end{array}$$

(9)

where $\mathbf{x}=[p,e_x,e_y,i_x,i_y]$ is a vector of modified equinoctial elements, and $\mathbf{g}(\mathbf{x},F,\mathbf{f})$ and $h(\mathbf{x},F,\mathbf{f})$ are the right-hand side of Equations (1)–(6).

In the current study, the optimal control problem between two points of phase space characterized by the elements and the eccentric longitudes $p_0,e_{x,0},e_{y,0},i_{x,0},i_{y,0},F_0$ and $p_1,e_{x,1},e_{y,1},i_{x,1},i_{y,1},F_1$ in a given time interval $[0,T]$ governed by Equations (8) is considered. In the current work, no other perturbations other than thrust are considered. The objective for minimization is $J={\int }_0^T{\left|\mathbf{f}\left(t\right)\right|}^{2}dt.$ This is the general minimum-energy transfer optimal control problem between two phase states in a given time [1]. It can be easily shown that in the case of an ideally regulated engine, this problem is equivalent to minimizing the fuel costs for a transfer between the two points in a given time.

Let us now proceed to the description of the averaging method used in this work.

3. The Averaging Technique and Approximate Solutions

3.1. Averaged Variational Equations in Modified Equinoctial Elements

In this section, we briefly describe the technique of obtaining the averaged variational equations and present the result of averaging. For that purpose, we first represent the perturbation acceleration by the Fourier series:

$$f_{r,c,n}={\alpha }_0^{r,c,n}+\sum _{k=1}^{\infty }({\alpha }_k^{r,c,n}coskF+{\beta }_k^{r,c,n}sinkF),$$

where ${\alpha }_k^{r,c,n}$ and ${\beta }_k^{r,c,n}$ are the Fourier series coefficients, and the upper script indicates one of the acceleration components. This representation allows us to average the right-hand sides of the variational equations on p, $e_x$, $e_y$, $i_x$, and $i_y$ over the mean longitude $\lambda =F+e_{x}cosF+e_{y}sinF$ over one revolution:

$$\begin{array}{cc}\hfill \frac{d\overline{\mathbf{x}}}{dt}& =\overline{\mathbf{g}}(\overline{\mathbf{x}},\mathbf{f})=\frac{1}{2\pi }{\int }_0^{2\pi }\mathbf{g}(\overline{\mathbf{x}},F\left(\lambda \right),\mathbf{f}\left(F\left(\lambda \right)\right))d\lambda =\hfill \\ & =\frac{1}{2\pi }{\int }_0^{2\pi }(1-{\overline{e}}_{x}cosF-{\overline{e}}_{y}sinF)\mathbf{g}(\overline{\mathbf{x}},F,\mathbf{f}\left(F\right))dF.\hfill \end{array}$$

(10)

Here, $\mathbf{x}$ is the vector of the modified equinoctial elements, and $\overline{\mathbf{x}}$ is the vector of their averaged values. Equations (1)–(5), averaged using the scheme (10), take the following form:

$$\begin{array}{cc}\hfill \frac{d\overline{p}}{dt}& =\frac{1}{{\varphi }^2}\sqrt{\frac{{\overline{p}}^3}{\mu }}\left\{(2+{\overline{e}}_x^2+{\overline{e}}_y^2){\alpha }_0^c-2{\overline{e}}_x{\alpha }_1^c-2{\overline{e}}_y{\beta }_1^c+\frac{{\overline{e}}_x^2-{\overline{e}}_y^2}{2}{\alpha }_2^c+{\overline{e}}_x{\overline{e}}_y{\beta }_2^c\right\},\hfill \\ \hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \frac{d{\overline{e}}_x}{dt}& =\sqrt{\frac{\overline{p}}{\mu }}\left\{-{\overline{e}}_y{\alpha }_0^r+{\overline{e}}_x{\overline{e}}_y\frac{1-\varphi }{2{\overline{e}}^2}{\alpha }_1^r+\frac{\varphi {\overline{e}}_x^2+{\overline{e}}_y^2}{2{\overline{e}}^2}{\beta }_1^r-\frac{3{\overline{e}}_x}{2}{\alpha }_0^c+\frac{{\overline{e}}_y^2(1+{\varphi }^2)+2\varphi {\overline{e}}_x^2}{2{\overline{e}}^2\varphi }{\alpha }_1^c-{\overline{e}}_x{\overline{e}}_y\frac{{(1-\varphi )}^2}{2{\overline{e}}^2\varphi }{\beta }_1^c-\right.\hfill \\ \hfill & -{\overline{e}}_x\frac{({\overline{e}}_x^2-{\overline{e}}_y^2)\varphi +2{\overline{e}}_y^2}{4{\overline{e}}^2\varphi }{\alpha }_2^c-{\overline{e}}_y\frac{2{\overline{e}}_x^2\varphi -{\overline{e}}_x^2+{\overline{e}}_y^2}{4{\overline{e}}^2\varphi }{\beta }_2^c+3{\overline{e}}_y\frac{{\overline{e}}_y{\overline{i}}_x-{\overline{e}}_x{\overline{i}}_y}{2{\varphi }^2}{\alpha }_0^n+\hfill \\ \hfill & +{\overline{e}}_y\frac{({\overline{e}}_x{\overline{i}}_y-{\overline{e}}_y{\overline{i}}_x)(2+\varphi ){\overline{e}}_x+\varphi (1+\varphi ){\overline{i}}_y}{2\left(1+\varphi \right){\varphi }^2}{\alpha }_1^n-{\overline{e}}_y\frac{({\overline{e}}_y{\overline{i}}_x-{\overline{e}}_x{\overline{i}}_y)(2+\varphi ){\overline{e}}_y+\varphi (1+\varphi ){\overline{i}}_x}{2\left(1+\varphi \right){\varphi }^2}{\beta }_1^n-\hfill \\ \hfill & \left.-{\overline{e}}_y\frac{({\overline{e}}_x{\overline{i}}_y+{\overline{e}}_y{\overline{i}}_x)\left(1+\varphi \right)-2{\overline{e}}_x{\overline{e}}_y({\overline{e}}_x{\overline{i}}_x+{\overline{e}}_y{\overline{i}}_y)}{4\left(1+\varphi \right){\varphi }^2}{\alpha }_2^n+{\overline{e}}_y\frac{\varphi (1+\varphi )({\overline{e}}_x{\overline{i}}_x-{\overline{e}}_y{\overline{i}}_y)+2{\overline{e}}_x{\overline{e}}_y({\overline{e}}_y{\overline{i}}_x-{\overline{e}}_x{\overline{i}}_y)}{4\left(1+\varphi \right){\varphi }^2}{\beta }_2^n\right\},\hfill \\ \hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \frac{d{\overline{e}}_y}{dt}& =\sqrt{\frac{\overline{p}}{\mu }}\left\{{\overline{e}}_x{\alpha }_0^r-\frac{\varphi {\overline{e}}_y^2+{\overline{e}}_x^2}{2{\overline{e}}^2}{\alpha }_1^r+{\overline{e}}_x{\overline{e}}_y\frac{\varphi -1}{2{\overline{e}}^2}{\beta }_1^r-\frac{3{\overline{e}}_y}{2}{\alpha }_0^c-{\overline{e}}_x{\overline{e}}_y\frac{{(1-\varphi )}^2}{2{\overline{e}}^2\varphi }{\alpha }_1^c+\frac{2\varphi {\overline{e}}_y^2+(1+{\varphi }^2){\overline{e}}_x^2}{2{\overline{e}}^2\varphi }{\beta }_1^c+\right.\hfill \\ \hfill & +{\overline{e}}_y\frac{{({\overline{e}}_y}^2-{\overline{e}}_x^2)\varphi +2{\overline{e}}_x^2}{4{\overline{e}}^2\varphi }{\alpha }_2^c-{\overline{e}}_x\frac{2{\overline{e}}_y^2\varphi +{\overline{e}}_x^2-{\overline{e}}_y^2}{4{\overline{e}}^2\varphi }{\beta }_2^c+3{\overline{e}}_x\frac{{\overline{e}}_x{\overline{i}}_y-{\overline{e}}_y{\overline{i}}_x}{2{\varphi }^2}{\alpha }_0^n-\hfill \\ \hfill & -{\overline{e}}_x\frac{({\overline{e}}_x{\overline{i}}_y-{\overline{e}}_y{\overline{i}}_x)\left(2+\varphi \right){\overline{e}}_x+\varphi (1+\varphi ){\overline{i}}_y}{2\left(1+\varphi \right){\varphi }^2}{\alpha }_1^n+{\overline{e}}_x\frac{({\overline{e}}_y{\overline{i}}_x-{\overline{e}}_x{\overline{i}}_y)\left(2+\varphi \right){\overline{e}}_y+\varphi (1+\varphi ){\overline{i}}_x}{2(1+\varphi ){\varphi }^2}{\beta }_1^n+\hfill \\ \hfill & \left.+{\overline{e}}_x\frac{({\overline{e}}_x{\overline{i}}_y+{\overline{e}}_y{\overline{i}}_x)(1+\varphi )-2{\overline{e}}_x{\overline{e}}_y({\overline{e}}_x{\overline{i}}_x+{\overline{e}}_y{\overline{i}}_y)}{4(1+\varphi ){\varphi }^2}{\alpha }_2^n-{\overline{e}}_x\frac{({\overline{e}}_x{\overline{i}}_x-{\overline{e}}_y{\overline{i}}_y)\varphi (1+\varphi )+2{\overline{e}}_x{\overline{e}}_y({\overline{e}}_y{\overline{i}}_x-{\overline{e}}_x{\overline{i}}_y)}{4(1+\varphi ){\varphi }^2}{\beta }_2^n\right\},\hfill \\ \hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \frac{d{\overline{i}}_x}{dt}& =\frac{1+{\overline{i}}_x^2+{\overline{i}}_y^2}{4}\sqrt{\frac{\overline{p}}{\mu }}\left\{-\frac{3{\overline{e}}_x}{2{\varphi }^2}{\alpha }_0^n+\frac{{\overline{e}}_x^2(2-{\varphi }^2)+\varphi {\overline{e}}_y^2}{{\overline{e}}^2{\varphi }^2}{\alpha }_1^n+{\overline{e}}_x{\overline{e}}_y\frac{2-{\varphi }^2-\varphi }{{\overline{e}}^2{\varphi }^2}{\beta }_1^n+{\overline{e}}_x\frac{{\overline{e}}_y^2-{\overline{e}}_x^2-2\varphi {\overline{e}}_y^2}{2{\overline{e}}^2{\varphi }^2}{\alpha }_2^n-\right.\hfill \\ \hfill & \left.-{\overline{e}}_y\frac{2{\overline{e}}_x^2-({\overline{e}}_x^2-{\overline{e}}_y^2)\varphi }{2{\overline{e}}^2{\varphi }^2}{\beta }_2^n\right\},\hfill \\ \hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{d{\overline{i}}_y}{dt}& =\frac{1+{\overline{i}}_x^2+{\overline{i}}_y^2}{4}\sqrt{\frac{\overline{p}}{\mu }}\left\{-\frac{3{\overline{e}}_y}{2{\varphi }^2}{\alpha }_0^n+{\overline{e}}_x{\overline{e}}_y\frac{2-\varphi -{\varphi }^2}{{\overline{e}}^2{\varphi }^2}{\alpha }_1^n+\frac{{\overline{e}}_y^2(2-{\varphi }^2)+{\overline{e}}_x^2\varphi }{{\overline{e}}^2{\varphi }^2}{\beta }_1^n-{\overline{e}}_y\frac{{\overline{e}}_x^2-{\overline{e}}_y^2-2{\overline{e}}_x^2\varphi }{2{\overline{e}}^2{\varphi }^2}{\alpha }_2^n-\right.\hfill \\ \hfill & \left.-{\overline{e}}_x\frac{2{\overline{e}}_y^2+({\overline{e}}_x^2-{\overline{e}}_y^2)\varphi }{2{\overline{e}}^2{\varphi }^2}{\beta }_2^n\right\},\hfill \end{array}$$

(15)

where $\overline{e}=\sqrt{{\overline{e}}_x^2+{\overline{e}}_y^2}$, $\overline{\varphi }=\sqrt{1-{\overline{e}}^2}$.

As one can see from the expressions above, the right-hand side of the equations depends only on the 13 Fourier coefficients: ${\alpha }_0^r$, ${\alpha }_1^r$, ${\beta }_1^r$, ${\alpha }_0^c$, ${\alpha }_1^c$, ${\beta }_1^c$, ${\alpha }_2^c$, ${\beta }_2^c$, ${\alpha }_0^n$, ${\alpha }_1^n$, ${\beta }_1^n$, ${\alpha }_2^n$, and ${\beta }_2^n$. Note that the major complexity in the right-hand sides of the equations is due to the terms proportional to ${\overline{e}}_x$ and ${\overline{e}}_y$. In the case of transfer between near-circular orbits, these terms can be neglected, and the averaged equations appear to be integrable. In the following subsection, the exact solutions to the simplified equations are given.

3.2. Closed-Form Averaged Solution of Orbital Motion Equations

In this section, the approximate closed-form solution of the averaged Equations (11)–(15) is presented for a low-thrust transfer between near-circular orbits.

A simple approximate relationship can be established between the error in the solutions of Equations (11)–(15) $|\Delta \overline{\mathbf{x}}|$ and the error in the calculation of their right-hand sides $|\Delta \overline{\mathbf{g}}|$:

$$|\Delta \overline{\mathbf{x}}|\approx |\Delta \overline{\mathbf{g}}|(t-t_0),$$

(16)

where $t_0$ is the initial moment of time. It can be noticed that in dimensionless units (with the unit of distance being the mean Earth radius $R_E=6371$ km and the unit of acceleration being the standard gravity $g_n=9.8067\mathrm{m}/{\mathrm{s}}^2$) for near-circular orbits ($e\le {10}^{-3}$) and low-thrust control ($\left|\mathbf{f}\right|\le {10}^{-4}$), the terms proportional to ${\overline{e}}_x$ and ${\overline{e}}_y$ in (11)–(15) are less than ${10}^{-6}$. Thus, neglecting them from equations yields an error of $|\Delta \overline{\mathbf{g}}|\le {10}^{-6}$ and, according to (16), over time intervals of several hundred revolutions ($t-t_0=2\pi ·160\approx {10}^3$), results in an error in solutions of the order ${10}^{-3}$ that can be considered acceptable.

After dropping the terms proportional to ${\overline{e}}_x$ and ${\overline{e}}_y$, Equations (11)–(15) have a simple and compact form:

$$\begin{array}{cc}\hfill \frac{d\overline{p}}{dt}& =2\sqrt{\frac{{\overline{p}}^3}{\mu }}{\alpha }_0^c,\hfill \\ \hfill \frac{d{\overline{e}}_x}{dt}& =\sqrt{\frac{\overline{p}}{\mu }}\left(\frac{{\beta }_1^r}{2}+{\alpha }_1^c\right),\hfill \\ \hfill \frac{d{\overline{e}}_y}{dt}& =\sqrt{\frac{\overline{p}}{\mu }}\left({\beta }_1^c-\frac{{\alpha }_1^r}{2}\right),\hfill \\ \hfill \frac{d{\overline{i}}_x}{dt}& =\sqrt{\frac{\overline{p}}{\mu }}\frac{1+{\overline{i}}_x^2+{\overline{i}}_y^2}{4}{\alpha }_1^n,\hfill \\ \hfill \frac{d{\overline{i}}_y}{dt}& =\sqrt{\frac{\overline{p}}{\mu }}\frac{1+{\overline{i}}_x^2+{\overline{i}}_y^2}{4}{\beta }_1^n.\hfill \end{array}$$

This system of equations is integrable; its solution, hereinafter called the zeroth-order solution, is presented below:

$$\begin{array}{cc}\hfill \overline{p}\left(\tau \right)& =\overline{p}\left(0\right)exp\left(2{\alpha }_0^c\tau \right),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\overline{e}}_x\left(\tau \right)& ={\overline{e}}_x\left(0\right)+\left(\frac{{\beta }_1^r}{2}+{\alpha }_1^c\right)\tau ,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\overline{e}}_y\left(\tau \right)& ={\overline{e}}_y\left(0\right)+\left({\beta }_1^c-\frac{{\alpha }_1^r}{2}\right)\tau ,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\overline{i}}_x\left(\tau \right)& =\frac{{\alpha }_1^n\rho tan\left(\gamma +\frac{\rho }{4}\tau \right)+{\beta }_1^{n}K}{{\rho }^2-K^2},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\overline{i}}_y\left(\tau \right)& =\frac{{\beta }_1^n\rho tan\left(\gamma +\frac{\rho }{4}\tau \right)-{\alpha }_1^{n}K}{{\rho }^2-K^2},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \tau \left(t\right)& =-\frac{1}{{\alpha }_0^c}ln\left(1-{\alpha }_0^c\sqrt{\frac{\overline{p}\left(t_0\right)}{\mu }}(t-t_0)\right),\hfill \end{array}$$

(22)

where

$$\begin{array}{cc}\hfill K& ={\beta }_1^n{\overline{i}}_x\left(0\right)-{\alpha }_1^n{\overline{i}}_y\left(0\right),\hfill \\ \hfill \rho & =\sqrt{{\left({\alpha }_1^n\right)}^2+{\left({\beta }_1^n\right)}^2+K^2},\hfill \\ \hfill \gamma & =arctan\left(\frac{{\alpha }_1^n{\overline{i}}_x\left(0\right)+{\beta }_1^n{\overline{i}}_y\left(0\right)}{\rho }\right).\hfill \end{array}$$

Thus, in the case of low eccentricity and perturbation magnitude, the evolution of the orbit can be described analytically, which reduces the time of trajectory propagation in comparison with the numerical integration of the initial averaged Equations (11)–(15). It is notable that the simplified equations depend on only seven Fourier coefficients: ${\alpha }_1^r$, ${\beta }_1^r$, ${\alpha }_0^c$, ${\alpha }_1^c$, ${\beta }_1^c$, ${\alpha }_1^n$, and ${\beta }_1^n$.

4. Multirevolution Low-Thrust Trajectory Optimization

The fact that the averaged Equations (11)–(15) are the functions of a finite number of Fourier coefficients allows one to parameterize the continuous perturbation acceleration and use direct approaches to solve an optimal control problem in averaged dynamics.

Let us utilize the spectral approach for low-thrust optimization problems. Assume that the problem of optimal control of the spacecraft’s orbital motion is defined as in Section 2. The solution to this problem is divided into two stages: first, solve the optimal control problem within the framework of the averaged dynamics with an averaged objective function; second, use the resulting solution as a starting approximation to solve the original problem in nonaveraged dynamics.

At the first stage, an optimal transfer between two orbits is sought. The orbits are defined each by five averaged modified equinoctial elements: $\overline{p}$, ${\overline{e}}_x$, ${\overline{e}}_y$, ${\overline{i}}_x$, and ${\overline{i}}_y$, which are equal to the given nonaveraged elements. The averaged orbital dynamics are defined by the expressions (11)–(15) or, in case of sufficiently low-thrust magnitudes and eccentricities, which is the case in this work, by Equations (17)–(22). The control is searched among periodic functions of F with seven nonzero Fourier coefficients, ${\alpha }_1^r$, ${\beta }_1^r$, ${\alpha }_0^c$, ${\alpha }_1^c$, ${\beta }_1^c$, ${\alpha }_1^n$, and ${\beta }_1^n$, and others equal to zero. Thus, we consider a class of control functions that directly affect the averaged dynamics of the spacecraft.

To obtain the averaged objective, first, let us average the integrand $f_r^2+f_c^2+f_n^2$:

$$\frac{1}{2\pi }{\int }_0^{2\pi }(f_r^2+f_c^2+f_n^2)d\lambda =\frac{1}{2\pi }{\int }_0^{2\pi }(1-{\overline{e}}_{x}cosF-{\overline{e}}_{y}sinF)(f_r^2+f_c^2+f_n^2)dF=$$

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=k_0+{\overline{e}}_x{k}_x+{\overline{e}}_y{k}_y,$$

where

$$k_0=\frac{1}{2}\left({\left|\alpha \right|}^2+{\left({\alpha }_0^c\right)}^2\right),k_x=-{\alpha }_0^c{\alpha }_1^c,k_y=-{\alpha }_0^c{\beta }_1^c.$$

Thus, the averaged objective is

$$\overline{J}=k_{0}T+k_x{\int }_0^T{\overline{e}}_{x}dt+k_y{\int }_0^T{\overline{e}}_{y}dt,$$

which can be obtained analytically using (17)–(22).

At the second stage, we revisit the initial optimal control problem and use the Fourier coefficients obtained from the averaged dynamics as our initial estimate. It is possible to add more terms to the Fourier series and assign zero values to the newly introduced higher-order coefficients.

5. Results

Let us demonstrate the technique on some examples. Two test problems are considered: the maneuvering in the vicinity of the geostationary orbit (GEO) and the orbit raising maneuver.

5.1. Maneuvering in the Vicinity of the Geostationary Orbit

First, we consider a transfer between two orbits close to the geostationary orbit (Case A). Such maneuvers can be performed to correct the orbit of geodetic and navigation satellites. The initial and target modified equinoctial elements are chosen as follows: $p\left(0\right)=42,500$ km, $e_x\left(0\right)=0.0007$, $e_y\left(0\right)=0.0009$, $i_x\left(0\right)=0.014$, $i_y\left(0\right)=0.022$, $p\left(T\right)=42,164$ km, $e_x\left(T\right)=0.0001$, $e_y\left(T\right)=0.0$, $i_x\left(T\right)=0.044$, and $i_y\left(T\right)=0.0$. The time of transfer $T=20$ days.

The eccentricity of both the initial and target orbits is close to zero; thus, the closed-form solutions (17)–(22) can be used to approximate spacecraft orbital dynamics. The control is parameterized by the seven Fourier coefficients: ${\alpha }_1^r$, ${\beta }_1^r$, ${\alpha }_0^c$, ${\alpha }_1^c$, ${\beta }_1^c$, ${\alpha }_1^n$, and ${\beta }_1^n$. The optimal control vector in both the averaging stage and the nonaveraging stage is sought by the sequential least squares programming (SLSQP) method. The values of the optimal thrust Fourier coefficients are obtained after 15 iterations for the first stage and 14 iterations for the second stage. Both sets of the optimal coefficients are presented in

Table 1. Optimal thrust Fourier coefficients in averaged and nonaveraged dynamics for Case A.

Coefficient	Averaged Model	Nonaveraged Model
	Radial, $\mathrm{mm}/{\mathrm{s}}^2$
${\alpha }_0^r$	0	−0.0001
${\alpha }_1^r$	0.0005	0.0005
${\beta }_1^r$	−0.0005	−0.0005
	Circumferential, $\mathrm{mm}/{\mathrm{s}}^2$
${\alpha }_0^c$	−0.0070	−0.0070
${\alpha }_1^c$	−0.0008	−0.0009
${\beta }_1^c$	−0.0013	−0.0013
${\alpha }_2^c$	0	0.0001
${\beta }_2^c$	0	−0.0002
	Normal, $\mathrm{mm}/{\mathrm{s}}^2$
${\alpha }_0^n$	0	−0.0001
${\alpha }_1^n$	0.2129	0.2136
${\beta }_1^n$	−0.1561	−0.1559
${\alpha }_2^n$	0	0.0008
${\beta }_2^n$	0	−0.0007

. Since some of the thrust coefficients are not present in (17)–(22), their optimal value is zero. The averaged and nonaveraged optimal thrust coefficients are placed on a single diagram in Figure 1.

During the maneuver, 20 revolutions are made around the central body. From Table 1 and Figure 1, it can be seen that the thrust coefficients not present in (17)–(22) are changed from zero to some nonzero optimal values. Their influence on the cost function is insignificant as their values are several orders of magnitude smaller than the largest coefficients.

The solution derived from the averaging stage is depicted in Figure 2, along with the nonaveraged solution extracted through numerical integration from the initial orbit using the discovered optimal control of the averaging stage. This figure indicates a minor discrepancy between the solutions in the averaged and nonaveraged versions of the problem, even in the absence of optimization in the nonaveraged problem formulation.

The elements of the final orbit, as determined by both the averaged and nonaveraged dynamics, are compared with the components of the target trajectory in

Table 2. Osculating modified equinoctial elements of the start, target, and final orbits for Case A.

	Start Orbit	Target Orbit	Final Orbit	Final Orbit
			at First Stage	at Second Stage
p, km	42,500	42,164	42,163.9	42,163.9993
$e_x$	$7·{10}^{-4}$	$1·{10}^{-4}$	$1.2·{10}^{-4}$	$9.94·{10}^{-5}$
$e_y$	$9·{10}^{-4}$	0	$-7·{10}^{-6}$	$2·{10}^{-7}$
$i_x$	$1.4·{10}^{-2}$	$4.4·{10}^{-2}$	$4.38·{10}^{-2}$	$4.3998·{10}^{-2}$
$i_y$	$2.2·{10}^{-2}$	0	$-2·{10}^{-5}$	$6·{10}^{-7}$

. The comparison reveals that the final elements obtained on the averaging stage are close to the target ones, and after the second stage, they become nearly identical to the target elements, showing convergence of the method.

5.2. Spiral Orbit Raising Maneuver

In the subsequent example (Case B), an optimization of a multirevolution spiral orbit-raising trajectory is carried out. The starting and target orbits are both equatorial circular orbits with $p\left(0\right)=20,000$ km and $p\left(T\right)=40,000$ km, respectively. The transfer time is set at $T=40$ days. Since the eccentricities of the initial and target orbits are zero, the spacecraft’s dynamics on the averaging stage are approximated using the zeroth-order solution. In this scenario, the optimal control vector can be calculated analytically and is solely circumferential (as indicated in

Table 3. Optimal thrust Fourier coefficients in averaged and nonaveraged dynamics for Case B.

Coefficient	Averaged Model	Nonaveraged Model
	Radial, $\mathrm{mm}/{\mathrm{s}}^2$
${\alpha }_0^r$	0	0.0002
${\alpha }_1^r$	0	0.0017
${\beta }_1^r$	0	0.0011
	Circumferential, $\mathrm{mm}/{\mathrm{s}}^2$
${\alpha }_0^c$	0.3783	0.3784
${\alpha }_1^c$	0	0.0019
${\beta }_1^c$	0	−0.0033
${\alpha }_2^c$	0	0.0003
${\beta }_2^c$	0	0.0001
	Normal, $\mathrm{mm}/{\mathrm{s}}^2$
${\alpha }_0^n$	0	−0.0001
${\alpha }_1^n$	0	0
${\beta }_1^n$	0	0
${\alpha }_2^n$	0	−0.0001
${\beta }_2^n$	0	−0.0001

). Using this solution as an initial estimate, the nonaveraged optimal thrust Fourier coefficients are obtained. Figure 3 provides a graphical comparison of the averaged and nonaveraged optimal thrust coefficients.

The number of revolutions of the optimal trajectory is 80. Table 3 and Figure 3 reveal that the optimal thrust coefficients undergo substantial changes due to the inaccuracies of the initial estimate. However, despite these modifications, the nonaveraged optimization process successfully converges at an optimal solution after 16 iterations of the SLSQP procedure.

Figure 4 illustrates the evolution of the averaged and nonaveraged osculating values of the modified equinoctial elements, influenced by the achieved optimal control of the averaging stage. The nonaveraged values of $e_x$ and $e_y$ experience substantial changes during the transfer, despite their corresponding averaged elements remaining at a value of zero. The final and target equinoctial elements are depicted in

Table 4. Osculating modified equinoctial elements of the start, target, final, and corrected final orbits for Case B.

	Start Orbit	Target Orbit	Final Orbit	Final Orbit
			at First Stage	at Second Stage
p, km	20,000	40,000	39,999	39,999.7
$e_x$	0.0	0.0	$-0.002$	$-3·{10}^{-6}$
$e_y$	0.0	0.0	$0.003$	$-7·{10}^{-6}$
$i_x$	0.0	0.0	$7·{10}^{-8}$	$-3·{10}^{-8}$
$i_y$	0.0	0.0	$2·{10}^{-8}$	$-8·{10}^{-8}$

. This example highlights that, in certain instances, the disparity between the averaged and nonaveraged trajectories can be substantial, necessitating corrections to the averaged optimal control in the nonaveraged scenario to meet the target orbit constraint within nonaveraged dynamics.

5.3. Comparison with Pontryagin Maximum Principle Control

It should be noted that the thrust accelerations obtained for Cases A and B are computed under the assumption that the control is a periodic function over the eccentric longitude. The question is whether the optimal control in a general class of piecewise-continuous function is close to a periodic function or not. To find it out, the obtained control is compared with the one calculated using the well-known Pontryagin maximum principle.

The optimal control in the class of piecewise-continuous functions is constructed as follows. The optimal control problem is posed in the classical Bolza formulation for a transfer between two fixed phase vectors for a fixed transfer time:

$$J=\frac{1}{2}{\int }_0^T{\left|\mathbf{f}\left(t\right)\right|}^{2}dt\to min,$$

$$\dot{\mathbf{r}}=\mathbf{v},\dot{\mathbf{v}}=-\mu \mathbf{r}/r^3+\mathbf{f},\mathbf{r}\left(0\right)={\mathbf{r}}_0,\mathbf{v}\left(0\right)={\mathbf{v}}_0,\mathbf{r}\left(T\right)={\mathbf{r}}_T,\mathbf{v}\left(T\right)={\mathbf{v}}_T,$$

where $\mathbf{r}$ and $\mathbf{v}$ are the position and velocity of the spacecraft; $\mathbf{f}$ is the thrust acceleration; T is a given time; and ${\mathbf{r}}_0$, ${\mathbf{v}}_0$, ${\mathbf{r}}_T$, and ${\mathbf{v}}_T$ are the initial and boundary constraints.

From a practical standpoint, the considered optimal control problem is equivalent to maximizing the spacecraft’s mass at the final time assuming a fixed jet propulsion power. In fact, since jet power equals $P=\left(fm\right)v_{\mathrm{ex}}/2$, where $v_{\mathrm{ex}}$ represents the exhaust velocity and $fm$ is the thrust force, then

$$\dot{m}=-\frac{fm}{v_{\mathrm{ex}}}=-\frac{f^2{m}^2}{2P}.$$

After integrating the derived equations, one gets

$$m\left(T\right)=\frac{m_0}{1+\frac{m_0}{P}J},$$

so minimizing J leads to maximizing $m\left(T\right)$.

The Hamiltonian function is written as follows:

$$H=-\frac{1}{2}{\left|\mathbf{f}\right|}^2+{\mathbf{p}}_r^T\mathbf{v}-{\mathbf{p}}_v^T\mathbf{r}/r^3+{\mathbf{p}}_v^T\mathbf{f},$$

where ${\mathbf{p}}_r$ and ${\mathbf{p}}_v$ are the conjugate variables. Maximizing H with respect to $\mathbf{f}$ gives $\mathbf{f}={\mathbf{p}}_v$, so the extended equations are

$$\dot{\mathbf{r}}=\mathbf{v},\dot{\mathbf{v}}=-\mu \mathbf{r}/r^3+{\mathbf{p}}_v,{\dot{\mathbf{p}}}_r={\mathbf{p}}_v/r^3-3\mathbf{r}{\mathbf{r}}^T{\mathbf{p}}_v/r^5,{\dot{\mathbf{p}}}_v={\mathbf{p}}_r.$$

The problem of finding an optimal solution is reduced to finding the initial values of the conjugate variables ${\mathbf{p}}_r\left(0\right)$ and ${\mathbf{p}}_v\left(0\right)$ such that propagating the above extended system of equations with $\mathbf{r}\left(0\right)={\mathbf{r}}_0$ and $\mathbf{v}\left(0\right)={\mathbf{v}}_0$ gives the target values $\mathbf{r}\left(T\right)={\mathbf{r}}_T$ and $\mathbf{v}\left(T\right)={\mathbf{v}}_T$.

The derived boundary value problem can be solved using the differential parameter continuation method with respect to the gravitational parameter, described, for example, in [32]. To utilize this method, it is necessary to determine the number of revolutions around the center of attraction. This is typically performed manually by managing the resulting fuel costs or maximum values of control variables. It is worth noting that the averaging method presented earlier does not require the specification of the number of revolutions. Moreover, while the differential continuation method is regular, it does involve integrating a system of differential equations, which includes the inversion of a matrix on the right-hand side. When dealing with a large number of revolutions, this method may face convergence issues, which are further complicated by the need to manually select the number of turns. This issue is not present in the proposed averaging method.

The Pontryagin maximum principle, in conjunction with the method of differential parameter continuation, was used to address the problems of Cases A and B, resulting in the acquisition of Pontryagin extremals (the control functions that meet the necessary conditions for optimality). The acceleration control profiles that were derived from the maximum principle and the two-staged optimization method are compared in Figure 5 and Figure 6 for Cases A and B, respectively. It is evident that in Case A, the accelerations achieved are almost the same. The resulting values of the objective are $J_1$ = 30,205 ${\mathrm{mm}}^2/{\mathrm{s}}^3$ for the presented technique and $J_2$ = 30,204 ${\mathrm{mm}}^2/{\mathrm{s}}^3$ for the maximum principle, the relative difference $(J_1-J_2)/J_1=0.003\%$. However, in Case B, while the accelerations appear to be similar, the spectrum of the Pontryagin optimal control is gradually changing over time. The objective values are also similar: 247,366 ${\mathrm{mm}}^2/{\mathrm{s}}^3$ for the two-stage optimization method and 247,365 ${\mathrm{mm}}^2/{\mathrm{s}}^3$ for the maximum principle with a relative difference equal to $0.0004\%$.

6. Conclusions and Further Research

This study presents a novel approach to generate an initial estimation for low-thrust optimization problems. This method utilizes the closed-form approximation of the averaged dynamics of a nearly circular orbit, enabling the conversion of the multirevolution orbital transfer problem into a moderately sized nonlinear programming problem. The solution to this problem is found to be near the optimal control. Applying this newly developed approach, two test energy-optimal control problems were successfully solved. It was determined that for transfers between closely situated near-circular orbits, the obtained control serves as a precise initial estimate: further optimization in the nonaveraged model only results in minor modifications.

Numerous aspects of the introduced method are being considered for additional investigation. For instance, the search for optimal low-thrust acceleration is conducted among functions with a constant Fourier spectrum. However, the energy-optimal thrust acceleration found in the nonaveraged model through indirect methods possesses a spectrum that progressively changes over time. In such scenarios, a potential modification to the introduced method could involve searching for thrust acceleration with a spectrum that is piecewise constant. This would cause the number of optimization variables to increase N-fold, but it would remain finite and relatively low.

Another question is the ability to account for various disturbances to the equations of motion: nonsphericity of the central body, other planets’ attraction, atmospheric drag, solar radiation pressure, etc. Accounting for them in the model of motion will remove the ability to analytically average the equations and thus parameterize the control. One of the ideas is to approximate disturbing accelerations with periodic functions like it is done for the low-thrust acceleration. This approximation might be quite accurate for disturbing accelerations with a slowly changing magnitude and direction, such as the gravitational attraction of other bodies, solar radiation pressure, and atmospheric drag. Unfortunately, nonsphericity gravity field effects that have a major effect on low orbits are hard to represent this way. The question of accounting for these disturbances is a subject of further research.