ORBIT TRANSFER METHOD FOR A SPACECRAFT USING A CONTINUOUS OR QUASI-CONTINUOUS THRUST AND EMBEDDED DRIVING SYSTEM FOR IMPLEMENTING SUCH A METHOD

Abstract

An orbit transfer method for a spacecraft using a continuous or quasi-continuous thrust propulsion, the method comprises: the acquisition, at least once in each half-revolution of the spacecraft, of measurements of its position and of its velocity; the computation of a thrust control function as a function of the measurements; and the driving of the thrust in accordance with the control law; wherein the control law is obtained from a Control-Lyapunov function using orbital parameters, preferably equinoctial, of the spacecraft, averaged over at least one half-revolution. An embedded driving system for a spacecraft for implementing such a method and a spacecraft equipped with the driving system are provided.

Claims

1. An orbit transfer method for a spacecraft using a continuous or quasi-continuous thrust propulsion, the method comprising: the acquisition, at least once in each revolution of the spacecraft, of measurements of its position and of its velocity; the computation of a thrust control function as a function of said measurements; and the driving of said thrust in accordance with said control law; wherein said control law is obtained from a control-Lyapunov function proportional to a term-weighted sum, each term being representative of a quadratic error between a measured orbital parameter of an orbit of the spacecraft averaged over at least one half-revolution and the corresponding orbital parameter of a final target orbit, normalized relative to the maximum value of the drift, averaged over at least one half-revolution, of said orbital parameter of said spacecraft.

2. The method according to claim 1, wherein said or each said orbital parameter is an equinoctial orbital parameter.

3. The method according to claim 1, wherein the acquisition of said position and velocity measurements is performed on board said spacecraft by means of a GNSS receiver.

4. The method according to claim 1, wherein said computation of a control function is performed by an embedded processor on board said spacecraft.

5. The method according to claim 1, wherein said control-Lyapunov function comprises at least one multiplying term consisting of a barrier function imposing a constraint of maximum or minimum altitude of said spacecraft.

6. The method according to claim 1, wherein the weights of said weighted sum are constant and non-negative.

7. The method according to claim 6, wherein the weights of said weighted sum are all equal.

8. The method according to claim 1, comprising a step of numerical optimization of the weights of said weighted sum, this step being performed by a computer on the ground before the start of the transfer and at least once during the transfer.

9. The method according to claim 1, also comprising, after an initial phase of the transfer: the computation of an estimation of a longitude encounter error on the final target orbit; the modification, in said control-Lyapunov function, of an orbital parameter of said final target orbit representing the half-major axis thereof or of a weighting coefficient of this parameter; and during a terminal phase of the transfer, the maintenance of a constant value of said orbital parameter and of weighting coefficient.

10. The method according to claim 1, wherein said propulsion is of electric type.

11. The method according to claim 1, wherein the driving of said thrust comprises the determination of at least its orientation.

12. An embedded system for driving a spacecraft comprising: a continuous or quasi-continuous thrust propulsion system; a GNSS receiver configured to acquire, at least once in each revolution of the spacecraft, measurements of its position and of its velocity; a processor programmed to compute a thrust control function as a function of said measurements and to drive said continuous or quasi-continuous thrust propulsion system in accordance with said control law; wherein said processor is programmed to compute said control law from a Control-Lyapunov function proportional to a term-weighted sum, each term being representative of a quadratic error between a measured orbital parameter of an orbit of the spacecraft averaged over at least one half-revolution and the corresponding orbital parameter of a final target orbit normalized relative to the maximum value of the rate of change of said orbital parameter, averaged over at least one half-revolution of said spacecraft.

13. The system according to claim 12, wherein said or each said orbital parameter is an equinoctial orbital parameter.

14. A spacecraft equipped with an embedded driving system according to claim 12.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Other features, details and advantages of the invention will emerge on reading the description given with reference to the attached drawings given by way of example and which represent, respectively:

[0043] FIG. 1, the orbit transfer of a satellite obtained by using a continuous or quasi-continuous thrust;

[0044] FIG. 2, the definition of the orbital parameters;

[0045] FIG. 3, the flow diagram of a method according to an embodiment of the invention;

[0046] FIG. 4, a simplified functional diagram of a spacecraft equipped with a driving system according to an embodiment of the invention;

[0047] FIGS. 5a to 5d, graphs illustrating the implementation of a driving method according to the prior art; and

[0048] FIGS. 6a to 6c and 7a, 7b, graphs illustrating the implementation of a driving method according to an embodiment of the invention.

DETAILED DESCRIPTION

[0049] As has been mentioned above, the invention uses a control law obtained from a control-Lyapunov function. It is therefore important to define such a function and its use to obtain a control law.

[0050] Consider a system governed by the differential equation

{dot over (x)}=f(x, u)   (1)

in which x is a state vector which, in the case considered here, represents the deviation between the measured orbital parameters of a spacecraft and the parameters of the target orbit, u is a control vector which defines the thrust (orientation and possibly intensity) and f(x,u) is the function, deriving from the laws of orbital mechanics, which expresses the time variation of the state vector as a function of the present value of this vector and of the control. A control-Lyapunov function V(x) is a continuously derivable function, strictly positive for any x except x=0, such as V(x=0)=0 and

∀x≠0 ∃u{dot over (V)}(x, u)=ΔV(x).f(x, u)<0   (2).

[0051] For the state vector to evolve towards its target x=0 (that is to say for the orbit of the spacecraft to tend towards the target orbit), V(x) must be minimized; it is therefore essential to choose a control u which renders its time derivative {dot over (V)}(x,u) as negative as possible. It is therefore natural to take

u=arg min.sub.{u}{dot over (V)}   (3).

[0052] The aim is therefore to find a control-Lyapunov function which provides, via the equation (3), a control law that is close to optimality. Optimality can be defined, for example, by a minimum transfer time, a minimal consumption of propellant or by a combination of these objectives.

[0053] Before proposing a form for the control-Lyapunov function it is a good idea to look into the choice of the state vector x. According to the prior art, this vector is defined from the five “conventional” orbital parameters, which are illustrated using FIG. 2: [0054] the half-major axis “a” of the orbit O (assumed circular or elliptical); [0055] the eccentricity “e” (ratio of the centre C-focus F distance to the half-major axis a; has the value 0 in the case of a circular orbit and is strictly between 0 and 1 for an elliptical orbit); [0056] the inclination “i” of the orbital plane PO relative to a reference plane PR, which can for example be the ecliptic or the equator; [0057] the longitude of the ascending node NA “Ω”, measured relative to a reference direction DR; [0058] the argument of the periastron “ω”, which is the angle formed by the line of the nodes NA-ND and the direction of the periastron DPA in the orbital plane.

[0059] These orbital parameters present the drawback of being ill-defined, and therefore of having singular movement equations (division by zero), for low eccentricities (e≈0) and for low inclinations (i≈0). For this reason, the invention uses so-called “equinoctial” orbital parameters, of which the movement equations are never singular and which are defined by:

[00001]$\begin{array}{cc}\left[\begin{array}{c}a\\ e_x\\ e_y\\ h_x\\ h_y\\ L\end{array}\right]=\left[\begin{array}{c}a\\ e.Math..Math.\cos \left(\Omega +\omega \right)\\ e.Math..Math.\sin \left(\Omega +\omega \right)\\ \tan \left(i/2\right).Math.\cos \left(\Omega \right)\\ \tan \left(i/2\right).Math.\sin \left(\Omega \right)\\ \Omega +\omega +v\end{array}\right]& \left(4\right)\end{array}$

in which v is the true anomaly, that is to say the angle between the direction of the periastron DPA and the line linking the centre C to the position of the spacecraft VS. It will be noted that the parameters e.sub.x, e.sub.y can be considered the components of an “eccentricity” vector of modulus “e” and having the direction of the perigee for polar angle, whereas the parameters h.sub.x, h.sub.y are the components of an “inclination” vector of modulus tan(i/2) and having the direction of the ascending node for polar angle.

[0060] The use of the equinoctial orbital parameters is not an essential feature of the invention. If the target orbit exhibits a not-inconsiderable eccentricity and inclination, it is also possible to use conventional orbital parameters.

[0061] The orbital parameters—conventional or equinoctial—are defined only for a Kepler orbit, which is not the case of a spacecraft subjected to a thrust. However, by knowing, at each instant, the position and the velocity of the spacecraft, it is possible to compute the parameters of its osculating orbit, that is to say the orbit the craft would follow if the propulsion was instantaneously cut and in the absence of any other disturbance. The position and the velocity of the spacecraft are generally known by virtue of the use of a GNSS navigation system—by direct measurement or interpolation between two successive measurements.

[0062] The invention does not use the orbital parameters measured (or, to be more precise, computed from the measurements) as such, but these parameters averaged over at least one half-period of revolution of the spacecraft. In effect, the measured orbital parameters oscillate at the orbital frequency, which is detrimental to the stability of the control; these oscillations are eliminated by the averaging operation.

[0063] A control-Lyapunov function according to an embodiment of the invention can be written as a weighted sum of average quadratic errors of the equinoctial orbital parameters. More specifically, it can be given by:

[00002]$\begin{array}{cc}V\left(\stackrel{\_}{a},{\stackrel{\_}{e}}_x,{\stackrel{\_}{e}}_y,{\stackrel{\_}{h}}_x,{\stackrel{\_}{h}}_y\right)={w_a\left(\frac{\stackrel{\_}{a}-{\stackrel{\_}{a}}_T}{{\stackrel{.}{\stackrel{\_}{a}}}_{\max }}\right)}^2+{w_{e_x}\left(\frac{{\stackrel{\_}{e}}_x-{\stackrel{\_}{e}}_{x,T}}{{\stackrel{.}{\stackrel{\_}{e}}}_{x,\max }}\right)}^2+{w_{e_y}\left(\frac{{\stackrel{\_}{e}}_y-{\stackrel{\_}{e}}_{y,T}}{{\stackrel{.}{\stackrel{\_}{e}}}_{y,\max }}\right)}^2+{w_{h_x}\left(\frac{{\stackrel{\_}{h}}_x-{\stackrel{\_}{h}}_{x,T}}{{\stackrel{.}{\stackrel{\_}{h}}}_{x,\max }}\right)}^2+{w_{h_y}\left(\frac{{\stackrel{\_}{h}}_y-{\stackrel{\_}{h}}_{y,T}}{{\stackrel{.}{\stackrel{\_}{h}}}_{y,\max }}\right)}^2& \left(5\right)\end{array}$

[0064] in which ā, ē.sub.x, ē.sub.y, h.sub.x, h.sub.y are the equinoctial orbital parameters of the orbit of the spacecraft averaged over at least one half-revolution, are the ā.sub.T, ē.sub.x,T, ē.sub.y, T.sub.T, h.sub.xT, h.sub.y,T are the equinoctial orbital parameters of the target orbit, {dot over (a)}.sub.max, ė.sub.x,max, ė.sub.y,max, {dot over (h)}.sub.x,max, {dot over (h)}.sub.y,max are the secular drifts (that is to say the average time derivatives) of the parameters ā, ē.sub.x, ē.sub.y, h.sub.x, h.sub.y obtained by application of the control law which maximizes these secular drifts or, in an equivalent manner, the orbital increment of each orbital parameter and w.sub.j(j=a, e.sub.x, e.sub.y, h.sub.x, h.sub.y) are weighting coefficients that are non-negative and preferably strictly positive. In the simplest embodiment, which nevertheless gives satisfactory results, these weighting coefficients can all be equal to one another, and notably taken to be equal to 1. In a variant, in order to get closer to the optimality control law, it is possible to proceed with a numerical optimization of these performance levels. This optimization can be performed on the ground, before the start of the transfer, for example by taking the transfer time as cost function to be minimized. It is also possible to repeat the optimization during the transfer with a very low repetition rate (for example once every three months), which entails transmitting the new optimized parameters to the embedded processor of the spacecraft.

[0065] In practice, to find the values {dot over (a)}.sub.max, ė.sub.x,max, ė.sub.y,max, {dot over (h)}.sub.x,max, {dot over (h)}.sub.y,max, the control law which maximizes the variation of each equinoctial orbital element for each orbital position L is integrated; then, each time derivative over the period considered and over all the values of L is averaged, leaving the other orbital elements constant. This integration can be done numerically or analytically.

[0066] The abovementioned article “Optimisation of low-thrust orbit transfers using the Q-Law for the initial guess” also recommends an optimization of the weighting coefficients of the Q-law but by using genetic algorithms which prove very cumbersome from the computational point of view. In the case of the invention, by virtue of the use of average orbital parameters, it is possible to use simpler nonlinear optimization techniques.

[0067] It is possible to impose altitude constraints (altitude of the periastron r.sub.p greater than or equal to a first threshold r.sub.p,min and/or altitude of the apastron r.sub.p less than or equal or a second threshold r.sub.a,max) by multiplying the function V(x) given by the equation (5) by corresponding barrier functions .sub.ε.sup.−, B.sub.68 .sup.+ characterized by a numerical smoothing parameter ε. A barrier function is a continuous (preferably derivable) function of which the value tends rapidly towards infinity in approaching a limit value, while remaining relatively flat far from this value. The constrained control-Lyapunov function can be written:

[00003]$\begin{array}{cc}{V}_{\mathrm{con}}\left(\stackrel{\_}{a},{\stackrel{\_}{e}}_x,{\stackrel{\_}{e}}_y,{\stackrel{\_}{h}}_x,{\stackrel{\_}{h}}_y\right)=B_{\mathrm{.Math.}}^{+}\left(\frac{r_p\left(\stackrel{\_}{a},{\stackrel{\_}{e}}_x,{\stackrel{\_}{e}}_y\right)}{r_{p,\min }}\right).Math.B_{\mathrm{.Math.}}^{-}\left(\frac{r_a\left(\stackrel{\_}{a},{\stackrel{\_}{e}}_x,{\stackrel{\_}{e}}_y\right)}{r_{a,\max }}\right).Math.V\left(\stackrel{\_}{a},{\stackrel{\_}{e}}_x,{\stackrel{\_}{e}}_y,{\stackrel{\_}{h}}_x,{\stackrel{\_}{h}}_y\right)& \left(6\right)\end{array}$

The use of a control-Lyapunov function expressed as a function of average orbital parameters offers several advantages:

[0068] the average parameters do not oscillate at the orbital frequency (unlike the non-averaged measured orbital parameters); that allows for a rapid integration of the differential equation (1), necessary for computing the control vector u via the equation (3) applied to the averaged dynamics system with a time step which can be several revolutions;

[0069] the desired trend of the thrust is smoother;

[0070] the time trend of the half-major axis is monotonic (see FIG. 6a), which makes it possible to implement a simple method to satisfy a terminal longitude encounter condition.

[0071] This method consists, on approaching the target orbit, in:

[0072] predicting the longitude of arrival on the target orbit, and estimating the encounter error relative to the desired longitude;

[0073] modifying, in the control-Lyapunov, function, the parameter ā.sub.T by dynamically computing a small deviation to the target to adjust the average drift in geographic longitude; and

[0074] before the end of the transfer, restoring the initial value of this parameter.

[0075] The idea is to provide the longitude of arrival and, if it does not correspond to the target, to modify the parameter ā.sub.T so as to correct the error. In an equivalent manner, that amounts to changing the numerical weight associated with the half-major axis in a ratio equal to {(ā-ā.sub.T.sup.1)/(ā-ā.sub.T)}.sup.2. This correction is dynamic in as much as it has to be recomputed several times along the transfer trajectory, because the prediction of the longitude of arrival is affected by an error which tends to decrease in time.

[0076] Generally, it is pointless to implement the correction more than one or two months before the predicted date of arrival on the target orbit, because the prediction errors would be too great. Furthermore, the correction must be stopped at least one or two weeks before the predicted date of arrival to avoid having the amplitude of the corrections diverge (when there is little time left, a significant modification ā.sub.T is necessary to even slightly change the longitude of arrival). In the final phase of the transfer, the last modified value of the parameter ā.sub.T or the last modified value of its weight is kept.

[0077] FIG. 3 shows a flow diagram of a method according to the invention.

[0078] The first step is to measure the position and the velocity of the spacecraft, performed by GNSS (or by telemetry, but the autonomy of the craft is then lost) at least once per revolution. That makes it possible to determine, at a plurality of instants, the (equinoctial) orbital parameters of the spacecraft, which are then averaged. The average parameters thus obtained are used to compute the control-Lyapunov function. The weighting coefficients of this function can be optimized periodically by a computer located on the ground, and transmitted to the embedded processor. Furthermore, a target parameter can be modified temporarily to correct an estimated longitude encounter error, this estimation being in turn computed from the RGNSS measurements. Next, the control-Lyapunov function is used to compute the control function u, which drives the continuous or quasi-continuous propulsion system.

[0079] In practice, a GNSS receiver acquires position and velocity measurements at a high rate (several times per minute), but a filtering is generally performed on these measurements to retain only a few acquisitions (typically between 1 and 4) per revolution.

[0080] FIG. 4 is a very simplified functional diagram of a spacecraft VS equipped with a driving system according to the invention. The driving system comprises: a GNSS receiver (reference GNSS) which supplies position and velocity measurements; optionally a receiver (reference RSS) which receives, from a ground station, updates of the weighting parameters of the control-Lyapunov function and/or updates of the parameters of the target orbit and/or of other controls; an embedded processor PE which receives the position and velocity signals from the GNSS receiver (and possibly the data originating from the RSS receiver) and which computes a control signal u; and a continuous or quasi-continuous thrust propulsion system (generally electrical) which receives and applies this control signal.

[0081] A method according to the invention has been tested for the case of the stationing of a geostationary satellite from an elliptical and inclined injection orbit. Table 1 below gives the (conventional) orbital parameters of the initial orbit and of the target orbit:

TABLE-US-00001 TABLE 1 Orbital parameters Initial orbit Target orbit Half-major axis (a) 24505.9 km 42165 km Eccentricity (e) 0.725 0.001 Inclination (i) 7.05 deg 0.05 deg Ascending node long (Ω) 0 deg free Argument of the perigee (ω) 0 deg free

[0082] The “free” parameters are processed by setting the corresponding weight w.sub.j to zero, or by setting an orbital target equal to the initial parameter (this second method works less well than the first if the natural disturbances in dynamics are considered).

[0083] The satellite has an initial weight of 2000 kg, uses a propellant (Xenon) of specific impulse of 2000 s and its electrical propulsion system has a thrust of 0.35 N. Minimal time solutions were considered, in which the thrust always takes its maximum value and only its orientation is driven.

[0084] Table 2 illustrates the performance levels obtained by using an “optimal” driving, in the sense that it minimizes the duration of the transfer, computed by means of the T_3D technique (see the article by T. Dargent cited above) and by the method of the invention (with unitary weighting coefficients). The performance metrics considered are the transfer duration (in days), the consumption of propellant (in kg) and the Delta-V.

TABLE-US-00002 Performance metric T_3D Invention Transfer duration 137.289 days 138.735 days Consumption of propellant 211.673 kg 213.903 kg Delta-V 2194.1 m/s 2218.5 m/s

[0085] These results are very satisfactory, because the method according to the invention brings about a very low cost overhead (1.05% for the transfer duration and the consumption of propellant, 1.01% for the Delta-V) compared to the optimal solution while being much less costly in terms of computation resources, which makes it possible to be implemented by an embedded processor. Also, these performance levels could be further improved by optimizing the weighting coefficients of the control-Lyapunov function.

[0086] The application of the Q-Law (see the abovementioned article “Techniques for designing many-revolution, electric-propulsion trajectories”, case B) leads to a significantly higher Xenon consumption: 221 kg. The use of an optimization by genetic algorithm makes it possible to reduce this consumption to 213 kg, but at the cost of a considerable increase in computational complexity.

[0087] FIGS. 5a to 5d make it possible to follow the history of the orbital transfer obtained by the T_3D method. More specifically: [0088] FIG. 5a shows the time trend of the half-major axis of the orbit; [0089] FIG. 5b shows the time trend of the eccentricity; [0090] FIG. 5c shows the time trend of the inclination; and [0091] FIG. 5d shows the time trend of the orbital radius which oscillates greatly with the half-orbital period. [0092] FIGS. 6a to 6c make it possible to follow the history of the orbital transfer obtained by the method according to the invention. More specifically: [0093] FIG. 6a shows the time trend of the radius at the apogee ra (highest dotted line curve), at the radius of the perigee rp (lowest dotted line curve) and of the half-major axis a (continuous line curve); it can be seen that the trend of the half-major axis is monotonic, which makes it possible to impose observance of a longitude encounter condition, as was explained above. The trend of the radius at the apogee, however, is not monotonic because it proves optimal in raising the apogee to make the correction of the inclination more effective. The unit of length used for the y axis, designated DU, corresponds to 10 000 km. [0094] FIG. 6b shows the time trend of the eccentricity; and [0095] FIG. 6c shows the time trend of the inclination.

[0096] The profiles of the half-major axis and of the eccentricity very closely resemble those of the optimal control solutions. On the other hand, the inclination profile is substantially different, particularly at the end of transfer. This means that there are numerous trajectories, mutually different but which are “quasi-optimal”. The method according to the invention makes it possible to find one of them.

[0097] FIG. 7a shows the trend of the control-Lyapunov function—in fact, of its square root, the dimension of which is that of a time, and which can be considered an approximation by excess of the remaining duration of the transfer. The dotted line represents the true remaining duration, determined a posteriori. FIG. 7b represents the derivatives of the control-Lyapunov function relative to the equinoctial orbital parameters, which make it possible to determine the direction of the thrust.