Use of the Lyapunov Functions for Calculating the Locally Optimal Control of a Thrust Vector during Low-Thrust Interorbital Transfer

This paper presents a method for locally optimal control of the thrust vector of the electric propulsion system (EPS) for a spacecraft that performs a multiturn interorbital transfer from the initial elliptical orbit into a geostationary orbit (GSO). The control represents the time dependences of the angles that characterize the EPS thrust vector orientation in space. Here, it is assumed that the EPS is always on. The proposed control algorithm belongs to the class of feedback control algorithms and is based on using the Lyapunov functions. Numerical examples are presented, which characterize the operability of the proposed control technique. Considerable attention is paid to the comparison of given solutions with the optimal solutions obtained within the framework of the maximum principle formalism.


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][2][3][4][5][6][7][8][9][10][11][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][16][17][18][19][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 (1) where is the spacecraft position vector, t is time, μ is the gravitational parameter of the Earth, P is the value of the EPS jet thrust. is the thrust function ( when the EPS is on, and when the EPS is off), m is the spacecraft mass, is the unit vector (ort) of the EPS thrust vector, w is the effective exhaust velocity of the EPS, is the matrix of transition from the ITRS coordinate system into the ICRS coordinate system, 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: where r E is the equatorial radius of the Earth, c n0 are coefficients at zonal harmonics, c nm and d nm are nonnormalized coefficients at tesseral and sectorial harmonics, is the mth derivative of the Legendre polynomial , , , ϕ 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 where is the gravitational parameter of the jth celestial body (index 1 denotes the Moon, and index 2 denotes the Sun), and is the radius-vector of the jth celestial body in the GECS.
To calculate celestial bodies' positions , 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: here, is the focal parameter; is eccentricity; is inclination; is the true anomaly; is the latitude argument; , , and 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 and , where is the jet acceleration magnitude, is the pitch angle, and is the yaw angle.
We will consider the motion without turning off the EPS; i.e., we will assume that . The interorbital transfer problem is formulated as follows. It is required to find functions and such that the spacecraft be transferred from the initial state into the final state In the case of transfer into the GSO that we are considering, = 42164 km, = 0, and = 0.
Next, we introduce into consideration the Lyapunov function in the following form: From the necessary and sufficient conditions for the minimum of function (4), which is twice differentiated with respect to and , we can obtain the law of controlling the jet acceleration vector, which pro-  ; L t k p k e k i vides, at each time moment, the minimum of the Lyapunov function's derivative: (5) here, 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: and , where 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: where and are the components of orts of the radial, transversal, and binormal in the GECS: and Interorbital transfer duration 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 , 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: (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: ; here, is the semimajor axis of the osculating transfer trajectory and is the true longitude of a spacecraft.
The following notations are introduced in system (7): We note that the equinoctial variables are related with classical elements in the following way: here, is the pericenter argument and is the longitude of the ascending node of an osculating trajectory.
The spacecraft with initial mass m 0 has to be transferred from the pericenter of some initial elliptical orbit, the ascending node longitude and the pericenter argument of which are zero, into the GSO for the minimum time.
The final conditions will take the form The thrust vector control should provide the minimum of an interorbital transfer time. Thus, we consider the problem of minimizing the functional: In this case, the Hamiltonian of the optimal control problem can be written as follows: (8) here, is the variable conjugate to the semimajor axis, is the variable conjugate to the true longitude, and (j = 1 … 4) are the variables conjugate to the equinoctial elements g 1 , …, g 4 .
The jet acceleration vector orientation is found from the condition of the maximum of Hamiltonian (8): (9) here, The conjugate variables included in the optimal control law (9)-(10) can be found from the system of the following form: (11) here, is the vector of phase variables and is the vector of conjugate variables.
It can be shown that, for the minimum-time action problem, we have 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 and time of interorbital transfer .
The point of entering the final orbit is not fixed; so, from the transversality condition, we have . The interorbital transfer time can be found from the condition . Thus, the vector of residuals at the right end of the transfer trajectory can be presented as follows: here, is the semimajor axis of the final orbit (GSO) and x is the vector of unknown parameters: Thus, the problem of searching for the optimal control is reduced to the numerical solution of the system of equations of the following form: To solve this system, we can use the Powell's hybrid method (the HYBRID algorithm [26]  ; 0 , 1  2  2  2  2  2  2  2  1  2  3  1  2  3   3  2  2  2  1  2  3 ; ; The operability of the proposed control technique, based on applying the Lyapunov functions, will be analyzed for several numerical examples of a lowthrust spacecraft transfer into the GSO. Let us compare the obtained results with the results found within the framework of solution of an optimal minimumtime 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.
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.
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.
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.
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.  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 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. and e i k k The most significant distinctions are observed in the evolutions of the ascending node longitude and pericenter argument (see Fig. 4).
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.
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.

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.

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/.