INTRODUCTION

Electric propulsion systems (EPSs) are widely used currently onboard spacecraft of various designation: from vehicles operating in low-Earth orbits to automatic interplanetary stations. Such propulsion systems are used both for trajectory correction and as cruise propulsion systems.

The high values of a thrust’s specific impulse possessed by EPSs allow significant savings of propellant (fuel) onboard the spacecraft. Because of this, interest in using such propulsion systems arose almost at the very beginning of the space age. Research in the field of analyzing the trajectories of EPS-equipped spacecraft has been carried out in many works beginning in the middle of the 20th century [1–12].

Most of these works were devoted to the problems of finding optimal (or quasi-optimal) control over an open loop. Such an approach to finding the control is quite legitimate for the tasks of designing spacecraft motion trajectories.

At the same time, to implement the control onboard the spacecraft, control algorithms should be developed, which, in accordance with the terminology of the automatic control theory, should be attributed to the class of feedback control algorithms. (By the way, such approaches can also be used for trajectory designing tasks.) One of the most effective approaches to the feedback control synthesis appears to have been presented in papers [13, 14].

The idea of using the Lyapunov functions to find the laws of thrust vector control for EPS-equipped spacecraft has also been rather fruitful, and many authors have used it [15–20].

The main advantage of the control algorithm used in this work is its simplicity. The operation of the proposed algorithm does not require large calculation resources. In addition, if necessary, it can be implemented onboard a spacecraft.

MATHEMATICAL MODEL OF SPACECRAFT MOTION

The motion of a spacecraft is considered in a nonrotating, geocentric equatorial coordinate system (GECS), the origin of which is located at the Earth’s center of masses and the axes of which are collinear to the axes of the International Celestial Reference System (ICRS). The International Terrestrial Reference System (ITRS) is used to calculate the disturbances from the noncentral nature of the Earth’s gravitational field.

The system of equations of the spacecraft motion in the GECS can be presented as

$$\begin{gathered} \frac{{{{d}^{2}}\boldsymbol{r}}}{{d{{t}^{2}}}} = - \frac{\mu }{{{{r}^{3}}}}\boldsymbol{r} + {\mathbf{Q}}\frac{{\partial R}}{{\partial {{r}_{g}}}} + {\mathbf{a}} + \frac{{\delta P}}{m}{{{\mathbf{e}}}_{p}}; \hfill \\ \frac{{dm}}{{dt}} = - \frac{{\delta P}}{w}; \hfill \\ \end{gathered} $$
(1)

where \(\boldsymbol{r}\) is the spacecraft position vector, t is time, μ is the gravitational parameter of the Earth, P is the value of the EPS jet thrust. \(\delta \) is the thrust function (\(\delta = 1\) when the EPS is on, and \(\delta = 0\) when the EPS is off), m is the spacecraft mass, \({{{\mathbf{e}}}_{p}}\) is the unit vector (ort) of the EPS thrust vector, w is the effective exhaust velocity of the EPS, \({\mathbf{Q}}\) is the matrix of transition from the ITRS coordinate system into the ICRS coordinate system, \({{\boldsymbol{r}}_{g}} = {{{\mathbf{Q}}}^{T}}r\) is the vector of spacecraft position in the ITRS coordinate system, R is the disturbing function caused by the noncentral nature of the Earth’s gravitational field, and a are disturbing accelerations from the attraction of the Moon and Sun.

The disturbing function caused by the noncentral nature of the Earth’s gravitational field is considered in the following form:

$$\begin{gathered} R = \frac{\mu }{r}\left[ {\sum\limits_{n = 2}^N {{{c}_{{n0}}}{{{\left( {\frac{{{{r}_{{\text{E}}}}}}{r}} \right)}}^{n}}P_{n}^{0}} } \right. \\ \left. { + \,\,\sum\limits_{n = 2}^M {{{{\left( {\frac{{{{r}_{{\text{E}}}}}}{r}} \right)}}^{n}}\sum\limits_{m = 1}^n {\left( {{{c}_{{nm}}}{{C}_{m}} + {{d}_{{nm}}}{{S}_{m}}} \right)P_{n}^{m}} } } \right], \\ \end{gathered} $$

where rE is the equatorial radius of the Earth, cn0 are coefficients at zonal harmonics, cnm and dnm are nonnormalized coefficients at tesseral and sectorial harmonics, \(P_{n}^{m} = {{d}^{m}}{{P}_{n}}{{\left( {\sin \varphi } \right)} \mathord{\left/ {\vphantom {{\left( {\sin \varphi } \right)} {d{{{\left( {\sin \varphi } \right)}}^{m}}}}} \right. \kern-0em} {d{{{\left( {\sin \varphi } \right)}}^{m}}}}\) is the mth derivative of the Legendre polynomial \({{P}_{n}}\left( {\sin \varphi } \right)\), \({{C}_{m}} = {{\cos }^{m}}\varphi \cos m\lambda \), \({{S}_{m}} = {{\cos }^{m}}\varphi \sin m\lambda \), φ is the geodetic latitude of a subsatellite point, λ is the longitude of a subsatellite point, and N and M are the order and degree of the used gravitational field model (the EGM96 model is used in the work, N = M = 70).

The disturbing accelerations from the attraction of the Moon and Sun are

$${\mathbf{a}} = \sum\limits_{j = 1}^2 {{{\mu }_{j}}\left( {\frac{{{{{\mathbf{r}}}_{j}} - {\mathbf{r}}}}{{{{{\left| {{{{\mathbf{r}}}_{j}} - {\mathbf{r}}} \right|}}^{3}}}} - \frac{{{{{\mathbf{r}}}_{j}}}}{{r_{j}^{{\mathbf{3}}}}}} \right)} ,$$

where \({{\mu }_{j}}\)is the gravitational parameter of the jth celestial body (index 1 denotes the Moon, and index 2 denotes the Sun), and \({{{\mathbf{r}}}_{j}}\) is the radius-vector of the jth celestial body in the GECS.

To calculate celestial bodies’ positions \({{{\mathbf{r}}}_{j}}\), the ephemeris software DE405 is used [21].

FEEDBACK CONTROL SYNTHESIS

The system of equations of the spacecraft disturbed motion in the osculating classical Keplerian elements is known [22]. In this case, the equations for the focal parameter of an orbit and for the eccentricity and inclination can be presented as follows:

$$\begin{gathered} \frac{{dp}}{{dt}} = 2p\sqrt {\frac{p}{\mu }} \frac{1}{{1 + e\cos \upsilon }}\left( {T + {{f}_{T}}} \right); \\ \frac{{de}}{{dt}} = \sqrt {\frac{p}{\mu }} \left( {\sin \upsilon \cdot \left( {S + {{f}_{S}}} \right)\frac{{^{{}}}}{{}}} \right. \\ \left. { + \,\,\frac{{e{{{\cos }}^{2}}\upsilon + 2\cos \upsilon + e}}{{1 + e\cos \upsilon }}\left( {T + {{f}_{T}}} \right)} \right); \\ \frac{{di}}{{dt}} = \sqrt {\frac{p}{\mu }} \frac{{\cos u}}{{1 + e\cos \upsilon }}\left( {W + {{f}_{W}}} \right); \\ \end{gathered} $$
(2)

here, \(p\) is the focal parameter; \(e\) is eccentricity; \(i\) is inclination; \(\upsilon \) is the true anomaly; \(u\) is the latitude argument; \({{f}_{S}}\),\({{f}_{T}}\), and \({{f}_{W}}\) are, respectively, the radial, transversal, and binormal projection of the jet acceleration acting on a spacecraft; and S, T, and W are, respectively, the radial, transversal, and binormal projection of the total disturbing acceleration from the other factors disturbing the trajectory (the noncentral nature of the Earth’s gravitational field, the Lunar–Solar disturbances, etc.).

Assuming that the EPS thrust vector is directed along the longitudinal axis of a spacecraft, we can write \({{f}_{S}} = f\sin \theta ;\) \({{f}_{T}} = f\cos \theta \cos \psi ;\) \({{f}_{W}} = f\cos \theta \sin \psi ,\) and \(f = {{\delta P} \mathord{\left/ {\vphantom {{\delta P} m}} \right. \kern-0em} m}\), where \(f = {{\delta P} \mathord{\left/ {\vphantom {{\delta P} m}} \right. \kern-0em} m}\) is the jet acceleration magnitude, \(\theta \) is the pitch angle, and \(\psi \) is the yaw angle.

We will consider the motion without turning off the EPS; i.e., we will assume that \(\delta \equiv 1\).

The interorbital transfer problem is formulated as follows. It is required to find functions \(\theta = \theta \left( t \right)\) and \(\psi = \psi \left( t \right)\) such that the spacecraft be transferred from the initial state \({{p}_{0}} = p\left( {{{t}_{0}}} \right);\) \({{e}_{0}} = e\left( {{{t}_{0}}} \right);\) \({{i}_{0}} = i\left( {{{t}_{0}}} \right)\) into the final state\({{p}_{f}} = p\left( {{{t}_{f}}} \right);\) \({{e}_{f}} = e\left( {{{t}_{f}}} \right);\) \({{i}_{f}} = i\left( {{{t}_{f}}} \right).\)

In the case of transfer into the GSO that we are considering, \({{p}_{f}}\)= 42 164 km, \({{e}_{f}}\)= 0, and \({{i}_{f}}\) = 0.

Next, we introduce into consideration the Lyapunov function in the following form:

$$L\left( t \right) = {{k}_{1}}\Delta {{p}^{2}} + {{k}_{2}}\Delta {{e}^{2}} + {{k}_{3}}\Delta {{i}^{2}};$$
(3)

where \({{k}_{1}},\;{{k}_{2}},\;{{k}_{2}}\) are some constant coefficients and \(\Delta p;\;\;\Delta e;\;\;\Delta i\) are the current dimensionless values of residuals in the focal parameter, eccentricity, and inclination (at the given point of the osculating transfer trajectory): \(\Delta p\left( t \right) = p\left( t \right) - {{p}_{f}},\) \(\Delta e\left( t \right) = e\left( t \right) - {{e}_{f}},\) and \(\Delta i\left( t \right) = i\left( t \right) - {{i}_{f}}.\)

Consider the time derivative of Lyapunov function (3):

$$\begin{gathered} g\left( {\theta ,\psi } \right) = \frac{{dL}}{{dt}} \\ = 2\left( {{{k}_{1}}\Delta p\frac{{d\Delta p}}{{dt}} + {{k}_{2}}\Delta e\frac{{d\Delta e}}{{dt}} + {{k}_{3}}\Delta i\frac{{d\Delta i}}{{dt}}} \right). \\ \end{gathered} $$
(4)

From the necessary and sufficient conditions for the minimum of function (4), which is twice differentiated with respect to \(\theta \) and \(\psi \), we can obtain the law of controlling the jet acceleration vector, which provides, at each time moment, the minimum of the Lyapunov function’s derivative:

$$\begin{gathered} {{f}_{S}} = f\frac{{ - {{\lambda }_{S}}}}{{\sqrt {\lambda _{T}^{2} + \lambda _{S}^{2} + \lambda _{W}^{2}} }};\;\,\,\,\,{{f}_{T}} = f\frac{{ - {{\lambda }_{T}}}}{{\sqrt {\lambda _{T}^{2} + \lambda _{S}^{2} + \lambda _{W}^{2}} }}; \\ {{f}_{W}} = f\frac{{ - {{\lambda }_{W}}}}{{\sqrt {\lambda _{T}^{2} + \lambda _{S}^{2} + \lambda _{W}^{2}} }}; \\ \end{gathered} $$
(5)

here,

$$\begin{gathered} {{\lambda }_{S}} = {{k}_{2}}\Delta e\sin \upsilon \left( {1 + \tilde {e}\cos \upsilon } \right); \\ {{\lambda }_{T}} = 2{{k}_{1}}\Delta p\tilde {p} + {{k}_{2}}\Delta e\left[ {\tilde {e}\left( {{{{\cos }}^{2}}\upsilon + 1} \right) + 2\cos \upsilon } \right]; \\ {{\lambda }_{W}} = {{k}_{3}}\Delta i\cos u;\,\,\,\tilde {p} = \Delta p + {{p}_{f}};\,\,\,\tilde {e} = \Delta e + {{e}_{f}}. \\ \end{gathered} $$

Assuming that the value of the focal parameter of the initial orbit differs from the required final value, and the initial orbit itself is elliptical and has a nonzero inclination, the coefficients, included in (3), can be introduced in the following manner: \({{k}_{1}} = {1 \mathord{\left/ {\vphantom {1 {{{{\left( {{{p}_{f}} - {{p}_{0}}} \right)}}^{2}}}}} \right. \kern-0em} {{{{\left( {{{p}_{f}} - {{p}_{0}}} \right)}}^{2}}}},\) \({{k}_{2}} = {{{{k}_{e}}} \mathord{\left/ {\vphantom {{{{k}_{e}}} {e_{0}^{2}}}} \right. \kern-0em} {e_{0}^{2}}},\) and \({{k}_{3}} = {{{{k}_{i}}} \mathord{\left/ {\vphantom {{{{k}_{i}}} {i_{0}^{2}}}} \right. \kern-0em} {i_{0}^{2}}}\) \(\Delta i\left( t \right) = i\left( t \right) - {{i}_{f}}\), where \({{k}_{e}},\;{{k}_{i}}\) are constant coefficients that characterize the weight of a residual with respect to the eccentricity and inclination, respectively.

It is not difficult to find the components of the thrust ort in the GECS:

$${{{\mathbf{e}}}_{p}} = \left( {\begin{array}{*{20}{c}} {{{r}_{x}}}&{{{n}_{x}}}&{{{b}_{x}}} \\ {{{r}_{y}}}&{{{n}_{y}}}&{{{b}_{y}}} \\ {{{r}_{z}}}&{{{n}_{z}}}&{{{b}_{z}}} \end{array}} \right) \times \left( {\begin{array}{*{20}{c}} {{{f}_{S}}} \\ {{{f}_{T}}} \\ {{{f}_{W}}} \end{array}} \right),$$

where \({{\left( {\begin{array}{*{20}{c}} {{{r}_{x}}}&{{{r}_{y}}}&{{{r}_{z}}} \end{array}} \right)}^{T}} = {\mathbf{r}_{0}};\) \({{\left( {\begin{array}{*{20}{c}} {{{n}_{x}}}&{{{n}_{y}}}&{{{n}_{z}}} \end{array}} \right)}^{T}} = {\mathbf{n}_{0}};\) and \({{\left( {\begin{array}{*{20}{c}} {{{b}_{x}}}&{{{b}_{y}}}&{{{b}_{z}}} \end{array}} \right)}^{T}} = {\mathbf{b}_{0}}\) are the components of orts of the radial, transversal, and binormal in the GECS: \({\mathbf{r}_{0}} = \frac{\mathbf{r}}{r},\) \({\mathbf{n}_{0}} = {\mathbf{b}_{0}} \times {\mathbf{r}_{0}},\) and \({\mathbf{b}_{0}} = \frac{{\mathbf{r} \times \dot {\mathbf{r}}}}{{\left| {\mathbf{r} \times \dot {\mathbf{r}}} \right|}}.\)

Interorbital transfer duration \({{t}_{f}}\) is determined by the time interval during which all parameters of the osculating transfer trajectory take the required values at the right end. At the same time, it should be noted that the orbital elements may come to their required values nonsimultaneously.

By varying weight coefficients \({{k}_{e}}\,\,{\text{and}}\,\,{{k}_{i}}\), it is possible to obtain different solutions to the interorbital transfer problem. All of them will differ in transfer duration, which, therefore, can be considered as some function of selected weight coefficients:

$${{t}_{f}} = {{t}_{f}}\left( {{{k}_{e}},{{k}_{i}}} \right).$$
(6)

By performing numerical minimization of (6) in the space of weight coefficients, it is possible to reach the decrease of the duration of a transport operation of the transfer and to ensure simultaneous fulfillment of all boundary conditions.

At the same time, it should be noted that the accuracy of fulfillment of boundary conditions at the right end of the trajectory, unfortunately, depends on the values of selected weight coefficients themselves. This results in the situation in which, in the process of numerical searching for the minimum of function (6), solutions appear that do not satisfy the required accuracy. These solutions need to be rejected. As a result, function (6) will have regions where it is not actually defined. Therefore, to find its minimum, it is highly desirable to use any nongradient searching techniques, such as those presented in papers [23, 24].

OPTIMAL CONTROL

The operability of the presented approach will be illustrated by several numerical examples. However, here, the natural question arises of how close the obtained solutions are to optimal ones, for example, in terms of the transfer duration criterion.

To answer this question and to find the optimal (in terms of the minimum-time action criterion) trajectories of transfer into the GSO, we have used the indirect method of optimizing the trajectories of motion of the EPS-equipped spacecraft. This method, which is based on applying the Pontryagin maximum principle formalism [25], allows us to reduce the problem of finding the optimal control to the boundary value problem for the system of ordinary differential equations.

At the first stage, we will analyze the motion of the EPS-equipped spacecraft in the central Newtonian field of the Earth’s gravity. In this case, it is convenient to write the equations of motion in the equinoctial elements that do not have singularities in the vicinity of zero eccentricity and inclination:

$$\begin{gathered} \frac{{dA}}{{dt}} = \frac{{2\sqrt {{{A}^{3}}} }}{{\sqrt {\gamma \mu } }}\left[ {\xi {{f}_{T}} + \left( {{{g}_{1}}\sin F - {{g}_{2}}\cos F} \right){{f}_{S}}} \right]; \\ \frac{{d{{g}_{1}}}}{{dt}} = \frac{1}{\xi }\sqrt {\frac{{A\gamma }}{\mu }} \\ \times \,\,\left\{ {\left[ {\left( {1 + \xi } \right)\cos F + {{g}_{1}}} \right]{{f}_{T}} + \xi \sin F{{f}_{S}} + {{g}_{2}}\eta {{f}_{W}}} \right\}; \\ \frac{{d{{g}_{2}}}}{{dt}} = \frac{1}{\xi }\sqrt {\frac{{A\gamma }}{\mu }} \\ \times \,\,\left\{ {\left[ {\left( {1 + \xi } \right)\sin F + {{g}_{2}}} \right]{{f}_{T}} - \xi \cos F{{f}_{S}} - {{g}_{1}}\eta {{f}_{W}}} \right\}; \\ \frac{{d{{g}_{3}}}}{{dt}} = - \sqrt {\frac{{A\gamma }}{\mu }} \frac{{\phi \cos F}}{{2\xi }}{{f}_{W}}; \\ \frac{{d{{g}_{4}}}}{{dt}} = - \sqrt {\frac{{A\gamma }}{\mu }} \frac{{\phi \sin F}}{{2\xi }}{{f}_{W}}; \\ \frac{{dF}}{{dt}} = \sqrt {\frac{\mu }{{{{{\left( {A\gamma } \right)}}^{3}}}}} {{\xi }^{2}} \\ - \,\,\frac{1}{\xi }\sqrt {\frac{{A\gamma }}{\mu }} \left( {{{g}_{3}}\sin F - {{g}_{4}}\cos F} \right){{f}_{W}}; \\ \frac{{dm}}{{dt}} = - \frac{{\delta P}}{w}; \\ \end{gathered} $$
(7)

here, \(A\) is the semimajor axis of the osculating transfer trajectory and \(F\) is the true longitude of a spacecraft.

The following notations are introduced in system (7):

$$\begin{gathered} \xi = 1 + {{g}_{1}}\cos F + {{g}_{2}}\sin F;\,\,\,\,\eta = {{g}_{3}}\sin F - {{g}_{4}}\cos F; \\ \phi = 1 + g_{3}^{2} + g_{4}^{2};\,\,\,\gamma = 1 - g_{1}^{2} - g_{2}^{2}. \\ \end{gathered} $$

We note that the equinoctial variables are related with classical elements in the following way:

$$\begin{gathered} {{g}_{1}} = e\cos \left( {\omega + \Omega } \right);\,\,\,\,{{g}_{2}} = e\sin \left( {\omega + \Omega } \right); \\ {{g}_{3}} = {\text{tan}}\frac{i}{2}\cos \Omega ;\,\,\,{{g}_{4}} = {\text{tan}}\frac{i}{2}\sin \Omega ; \\ F = \Omega + \omega + \upsilon ; \\ \end{gathered} $$

here, \(\omega \) is the pericenter argument and \(\Omega \) is the longitude of the ascending node of an osculating trajectory.

The spacecraft with initial mass m0 has to be transferred from the pericenter of some initial elliptical orbit, the ascending node longitude and the pericenter argument of which are zero,

$$\begin{gathered} A\left( {{{t}_{0}}} \right) = {{A}_{0}};\,\,\,\,{{g}_{1}}\left( {{{t}_{0}}} \right) = {{g}_{{10}}};\,\,\,\,{{g}_{2}}\left( {{{t}_{0}}} \right) = {{g}_{{20}}}, \\ {{g}_{3}}\left( {{{t}_{0}}} \right) = {{g}_{{30}}};\,\,\,{{g}_{4}}\left( {{{t}_{0}}} \right) = {{g}_{{40}}};\,\,\,F\left( {{{t}_{0}}} \right) = 0, \\ \end{gathered} $$

into the GSO for the minimum time.

The final conditions will take the form

$$\begin{gathered} A\left( {{{t}_{f}}} \right) = 42\,164\,\,{\text{km}};\,\,\,\,{{g}_{1}}\left( {{{t}_{f}}} \right) = 0; \\ {{g}_{2}}\left( {{{t}_{f}}} \right) = 0;\,\,\,\,{{g}_{3}}\left( {{{t}_{f}}} \right) = 0;\,\,\,{{g}_{4}}\left( {{{t}_{f}}} \right) = 0. \\ \end{gathered} $$

The thrust vector control should provide the minimum of an interorbital transfer time. Thus, we consider the problem of minimizing the functional: \(J = \int_{{{t}_{0}}}^{{{t}_{f}}} {dt \to \min .} \)

In this case, the Hamiltonian of the optimal control problem can be written as follows:

$$H = - 1 + {{\lambda }_{A}}\frac{{dA}}{{dt}} + {{\lambda }_{F}}\frac{{dF}}{{dt}} + \sum\limits_{j = 1}^4 {{{\lambda }_{j}}} \frac{{d{{g}_{j}}}}{{dt}};$$
(8)

here, \({{\lambda }_{A}}\) is the variable conjugate to the semimajor axis, \({{\lambda }_{F}}\) is the variable conjugate to the true longitude, and \({{\lambda }_{j}}\) (j = 1 … 4) are the variables conjugate to the equinoctial elements g1, …, g4.

The jet acceleration vector orientation is found from the condition of the maximum of Hamiltonian (8):

$$\begin{gathered} {{f}_{T}} = f\frac{{{{b}_{1}}}}{{\sqrt {b_{1}^{2} + b_{2}^{2} + b_{3}^{2}} }};\,\,\,{{f}_{S}} = f\frac{{{{b}_{2}}}}{{\sqrt {b_{1}^{2} + b_{2}^{2} + b_{3}^{2}} }}; \\ {{f}_{W}} = f\frac{{{{b}_{3}}}}{{\sqrt {b_{1}^{2} + b_{2}^{2} + b_{3}^{2}} }}; \\ \end{gathered} $$
(9)

here,

$$\begin{gathered} {{b}_{1}} = \frac{{2{{\lambda }_{A}}A\xi }}{\gamma } + {{\lambda }_{1}}\left( {\cos F + \frac{{{{g}_{1}} + \cos F}}{\xi }} \right) \\ + \,\,{{\lambda }_{2}}\left( {\sin F + \frac{{{{g}_{2}} + \sin F}}{\xi }} \right). \\ {{b}_{2}} = \frac{{2{{\lambda }_{A}}A}}{\gamma }\left( {{{g}_{1}}\sin F + {{g}_{2}}\cos F} \right) \\ + \,\,{{\lambda }_{1}}\sin F - {{\lambda }_{2}}\cos F. \\ {{b}_{3}} = \frac{1}{\gamma }\left[ {\eta \left( {{{\lambda }_{1}}{{g}_{2}} - {{\lambda }_{2}}{{g}_{1}}} \right) - \frac{\phi }{2}\left( {{{\lambda }_{3}}\cos F + {{\lambda }_{4}}\sin F} \right)} \right]. \\ \end{gathered} $$
(10)

The conjugate variables included in the optimal control law (9)–(10) can be found from the system of the following form:

$$\frac{{d \boldsymbol{\lambda} }}{{dt}} = - \frac{{\partial H}}{{\partial x}};$$
(11)

here, \(\boldsymbol{x} = {{\left( {A\;{{g}_{1}}\;{{g}_{2}}\;{{g}_{3}}\;{{g}_{4}}\;F} \right)}^{T}}\) is the vector of phase variables and \(\boldsymbol{\lambda} = {{\left( {{{\lambda }_{A}}\;{{\lambda }_{1}}\;{{\lambda }_{2}}\;{{\lambda }_{3}}\;{{\lambda }_{4}}\;{{\lambda }_{F}}} \right)}^{T}}\)is the vector of conjugate variables.

It can be shown that, for the minimum-time action problem, we have \(\delta = 1\) on the entire trajectory; that is, the spacecraft moves with the EPS turned on permanently.

Systems (7) and (11), when substituting control law (9) into (7), jointly form the system of equations of optimal motion of a spacecraft. To integrate this system, it is necessary to find the vector of conjugate variables at initial time moment \({{\boldsymbol{\lambda} }_{0}} = \boldsymbol{\lambda} \left( {{{t}_{0}}} \right)\) and time of interorbital transfer \({{t}_{f}}\).

The point of entering the final orbit is not fixed; so, from the transversality condition, we have \({{\lambda }_{F}}\left( {{{t}_{f}}} \right) = 0\). The interorbital transfer time can be found from the condition \(H\left( {{{t}_{f}}} \right) = 0\). Thus, the vector of residuals at the right end of the transfer trajectory can be presented as follows:

$$\boldsymbol{g}(\boldsymbol{x}) = {{\left[ {A\left( {{{t}_{f}}} \right) - A^{*}{{g}_{1}}\left( {{{t}_{f}}} \right){{g}_{2}}\left( {{{t}_{f}}} \right){{g}_{3}}\left( {{{t}_{f}}} \right){{g}_{4}}\left( {{{t}_{f}}} \right)\;{{\lambda }_{F}}\left( {{{t}_{f}}} \right)H\left( {{{t}_{f}}} \right)} \right]}^{T}};$$

here, \(A{\text{*}}\) is the semimajor axis of the final orbit (GSO) and x is the vector of unknown parameters:

$$\boldsymbol{x} = \left( \begin{gathered} \boldsymbol{\lambda} \left( {{{t}_{0}}} \right) \hfill \\ {{t}_{f}} \hfill \\ \end{gathered} \right).$$

Thus, the problem of searching for the optimal control is reduced to the numerical solution of the system of equations of the following form: \(\boldsymbol{g}(\boldsymbol{x}) = 0.\)

To solve this system, we can use the Powell’s hybrid method (the HYBRID algorithm [26]).

NUMERICAL RESULTS

The operability of the proposed control technique, based on applying the Lyapunov functions, will be analyzed for several numerical examples of a low-thrust spacecraft transfer into the GSO. Let us compare the obtained results with the results found within the framework of solution of an optimal minimum-time action problem. For both control versions, the spacecraft motion was analyzed within the framework of the model of the central Newtonian field of the Earth’s gravity.

As an example, we consider the transport system based on the Soyuz-2-1B launch vehicle (LV) and the Fregat upper stage (US). The scheme of payload insertion into the GSO is as follows. The LV provides the head unit insertion into the parking circular orbit. The head unit includes the US, the payload adapter, and the payload itself, which represents the spacecraft intended for inserting into the GSO. With the help of the US, the head unit is transferred into some intermediate elliptical orbit, the pericenter argument of which is assumed to be zero. On this orbit, the US is separated from the spacecraft together with the payload adapter. Further insertion of a spacecraft is accomplished under an effect of the EPS thrust force.

The mass of the head unit, inserted into the parking circular orbit with altitude of 200 km and inclination of 51.7° is 8320 kg when launching from the Vostochny (Eastern) cosmodrome. (https://www.samspace.ru/ products/launch_vehicles/rn_soyuz_2/, with the date of access being June 8, 2020). The final mass of the US is 1050 kg, and the specific impulse of thrust of its cruise propulsion system is 333.2 s [27]. The payload adapter mass is assumed to be 50 kg. The spacecraft’s EPS includes two SPD-100D Hall’s stationary plasma thrusters, which operate simultaneously. The thrust of one engine is 90 mN, and the specific impulse of thrust is 1740 s (https://fakel-russia.com/index. php/ru/produktsiya, with the date of access being June 8, 2020).

The versions of intermediate elliptical orbits that provide various times of transfer into the GSO are presented in Table 1. The mass of the spacecraft delivered into the intermediate orbit by the US was estimated under the assumption that the interorbital transfer is a two-pulse apsidal one. The value of gravitational losses and losses for control at implementing the first impulse is 2.5%. The second impulse is implemented perfectly.

Table 1. Versions of intermediate elliptical orbits
Full size table

The EPS thrust vector control for the interorbital transfer from the indicated intermediate orbits into the GSO was found by applying two approaches described above: the feedback control based on the Lyapunov functions and optimal control within the framework of Pontryagin’s maximum principle.

Table 2 presents the results of solution of interorbital transfer problems within the framework of the first approach. The numbers of solutions in this table correspond to the numbers of intermediate orbits from Table 1. The exit from the integration of spacecraft’s motion equations took place when the required value of the semimajor axis was reached. The residuals in the eccentricity and inclination implemented at the final time moment are also presented in Table 2.

Table 2. Results of solution of the problems of interorbital transfer into the GSO using feedback control based on the Lyapunov functions
Full size table

Now we estimate how close the obtained solutions are to optimal ones. Table 3 presents the main results of calculating the section of transfer from the intermediate orbit into the GSO for the two approaches under consideration. The numbers of solutions in the table correspond to the numbers of intermediate orbits from Table 1.

Table 3. Main results of calculating the section of transfer from the intermediate orbit into the GSO for two approaches under consideration
Full size table

It can be seen from the analysis of Table 3 that the solutions obtained within the framework of the proposed approach were very close to optimal in terms of the transfer duration and, accordingly, in terms of the final mass of a spacecraft. Indeed, the transfer duration in the optimal solutions is only 2–4% less than in the solutions obtained with the use of the discussed control technique based on the Lyapunov functions. The final mass of the spacecraft delivered into the GSO is almost the same for both control versions: the distinction does not exceed 4.51 kg.

We will also analyze the evolution of basic orbital elements of the transfer trajectory for the case of using feedback control based on the Lyapunov functions and for the case of optimal control. As an example, we consider the fourth solution (in Tables 1–3, it is highlighted in italics). The duration of transfer into the GSO is about 6 months in this case.

Figure 1 shows the time dependences of a semimajor axis, eccentricity, and inclination of the osculating transfer trajectory for the optimal control and for the case of using the feedback control on the basis of the Lyapunov functions.

Fig. 1.
figure 1

Time dependences of the transfer trajectory elements.

Full size image

As we can see, the character of changing the eccentricity and inclination is, generally, similar in both solutions, with the strongest difference being observed in the character of changing the semimajor axis of an orbit.

Now, we analyze the resulting control for both versions of solutions. Figures 2 and 3 present the dependences of the angles of declination and right ascension of the EPS thrust vector on the flight time in the GECS.

Fig. 2.
figure 2

Dependences of the angle of declination of the EPS thrust vector on the transfer time.

Full size image
Fig. 3.
figure 3

Dependences of the angles of right ascension of the EPS thrust vector on the transfer time.

Full size image

The changes of declination and right ascension angles have an oscillatory character in both types of solutions. It can be noted that the character of control in the optimal solution differs from that within the framework of the proposed control approach. Despite this, the transfer durations, both in the optimal solutions and in the solutions corresponding to feedback control based on applying the Lyapunov function are quite close in magnitude.

DISTURBED MOTION OF A SPACECRAFT

One advantage of the proposed approach to the EPS thrust vector control is its simplicity of adaptation to the case of disturbed motion. In the final section of the paper, we will give an example of implementing the control technique under consideration within the framework of the updated model of motion (1) taking into account the gravitational effect on the trajectory of transfer from the Moon and Sun, as well as the noncentral nature of the Earth’s gravitational field.

We now compare the solutions obtained within the framework of the model of the central field of the Earth’s gravity and within the framework of the updated model. As an example, we will consider the above-analyzed version of the interorbital transfer into the GSO (version 4 from Tables 1–3).

In the case of disturbed motion, it is necessary to specify the initial time moment to calculate the disturbing accelerations. Such a moment, i.e., the time moment of passage through the pericenter of intermediate orbit 4, is considered to be 00.00.00 UTC on January 1, 2025 (see Table 1). Weight coefficients \({{k}_{e}}\,\,{\text{and}}\,\,{{k}_{i}}\) included in the control law are the same as for the case of undisturbed motion (they are presented in Table 2).

Analysis of disturbed motion of a spacecraft has shown that the time dependences of change of a semimajor axis, eccentricity and inclination, during the transfer, do not virtually differ from similar dependences for the case of motion in the central field.

The most significant distinctions are observed in the evolutions of the ascending node longitude and pericenter argument (see Fig. 4).

Fig. 4.
figure 4

Evolution of the longitude of the ascending node.

Full size image

The character of control in the disturbed and undisturbed solution also turns out to be close; however, in the case of disturbed motion, one can observe the shift in the phase of oscillations of the angles of declination and right ascension of the EPS thrust vector. This phase shift gradually grows towards the end of the flight. As an example, Fig. 5 shows the interval of 150–155 days of flight.

Fig. 5.
figure 5

Dependences of the angles of declination and the EPS thrust vector on the transfer time for the cases of disturbed and undisturbed motion.

Full size image

The inclusion of disturbances into the model led to a slight increase of the transfer duration (up to 186.82 days vs. 185.43 days of transfer in the central field). This resulted in a certain decrease of the final mass of a spacecraft, which, however, turned out to be insignificant. Figure 6 shows the projections of the obtained transfer trajectory in the GECS.

Fig. 6.
figure 6

Transfer trajectory projections in the GECS.

Full size image

CONCLUSIONS

So, within the framework of this work, a technique has been proposed for controlling the thrust vector of a spacecraft performing interorbital transfer into the GSO with the use of an EPS. This control technique is rather simple and is based on using the Lyapunov functions. Within the framework of the considered class of interorbital transfer problems, this technique makes it possible to find the solutions that are very close to optimal in terms of the transfer duration.

Another advantage of the proposed approach is the fact that it can easily be used for accurate models of spacecraft motion with a various set of disturbing factors.

Due to the fact that the proposed control law is specified as a function of the current phase state of a spacecraft, this control belongs to the class of feedback control algorithms, while the algorithm itself does not require large calculational resources and can be applied onboard the spacecraft.