Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.kiam1.rssi.ru/pubs/AAS_09-103.pdf Äàòà èçìåíåíèÿ: Fri Sep 5 11:21:20 2014 Äàòà èíäåêñèðîâàíèÿ: Sat Apr 9 22:26:01 2016 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï

AAS 09-103

ALGORITHM OF AUTOMATIC DETECTION AND ANA LYSIS OF NON-EVOLUTIONA RY CHANGES IN ORBI TAL MOTION OF GEOCENTRIC OBJECTS
Sergey Kamensky*, Andrey Tuchin , Victor Stepanyants, Kyle T. Alfriend§
The goal of this work is develop ing the methods and algor ithms for automatic detection and analysis of non-evolu tionary changes in orbital motion of geocentric objects caused by man euvers and o ther events (but not by natural perturbations). The task is formu lated for three k inds of non-ev olution ary chang es, namely th e sing le impu lse man euver, two impulses maneu ver and small continuous thrust. The methods and algorithms ar e developed. Examp les of the algorithm work are given .

INTRODUCTI ON Automatic determination and adequate characterization of the man euvers performed by spacecraft in near-Earth orbits is a significant compon ent of the general p roblem o f tracking Earth satellites and catalog maintenance. This paper presents the techniques and algorithms for the detection and analysis of non-evolutionary changes of orbital motion of satellites resulting from maneuvers. The p roblem is formulated for three types of n on-evolutionary change of orbital parameters: one burn man euvers, t wo burn man euvers and lo w thrust man euvers. Identification of the non-evolutionary change of the orbit requires rather good knowledge of the satellite orbital motion, which can be acquired by processing of the measurements obtained by the sensors. This paper assumes that the satellite orbits are already determined before the maneu ver and after. In case the measurements were continuous we could consider only one burn and low thrust maneuvers. The techniques described in this work use a nu mber of classical problems . The paper describes the techniques for solving them wi th adjustments required by the go al task. The problem of the determination of the orbital transfers for given initial and final orbits has an infinite number of solutions. Selection of the most probable scheme of maneuvering is based on the evaluation of the characteristic velocity. Selecti on of one of three (one burn, two burn or low thrust) v ariants of man euvering is bas ed on the following criteria: We assume a one burn transfer wh en the orbits before and after the maneuver intersect and the value of the characteristic velocity does not ex ceed a sp ecified limit.

* §

Design man ager, "Vympel" Corporation, Mo scow , Russia Leading research er, "Vymp el" Corporation, Moscow, Russia Leading research er, "Vymp el" Corporation, Moscow, Russia TEES Distinguished Research Ch air Pro fesso r, Dept. of Aerosp ace Engineering, T exas A&M University , USA

1


The t wo burn of the minimum performed using Lambert problem

transfer is al ways possible. The times of the burns are selected by the condition of the characteristic velocity. The ev aluation of the two burn maneuver can be two techniques: by using the matrices of the partial derivatives or by using the depending on the ch aract er of orbital motion.

The lo w thrust maneuv er is evaluat ed in case we kno w the spacecraft is equipped with a lo w thrust engine. The first section of the report considers the mathematical basis of the algorithms for characterization of the maneuvers. We consider the classical problem of orbit determination given two positions and the time for the transfer bet ween them ­ the Lambert's problem and the problem of evaluation of the required characteristic velocity for coplanar transfers. Th e algorithm for the Lambert's problem is based on the techniques suggested by the Russian scientist Subbotin1,2,3. The geo metric method is used for evaluation of energy supply required for coplanar maneuvers. The second section of the report considers the evaluation of the impulse of the burn when we know the orbital paramet ers before and after the burn. In the case of a one burn maneuver the tracks before and after the burn must intersect in the point where the burn was applied. Thus the analysis of the cause of the change of the orbital motion should include determination of the point of minimu m distance bet ween the tracks before and aft er the burn . The resulting vector of residuals between orbital positions before and after the burn should be compared with the error of position prediction based on orbital parameters before and after the burn. The third section of the repo rt considers the case wh en the parameters of t wo orbits are kno wn and two burns perform the transfer between them. The algorithm for determination of the times of performing the burns and ev aluation of the respective characteristic velocities is described. The fou rth section describes the algorithm for determi nation the time interval for the work of low thrust engine and the generated acceleration provid ing the required orbital transfer.

MATHEMATI CAL BACKGROUND Dev elopment of the algorithms for evaluation of o rbital man euvers requires s everal classical methods. Description of these methods is given by Vallado 4. In particular this book includes a description of the methods of orbit determination for t wo given positions and the time interval for orbital transfer (Lamb ert Problem) and the analysis of coplanar maneuvers . This section of the paper considers these cl assical problems . The algorithm for Lambert's problem is based on the techniques suggested by Subbotin1,2,3. The geo metric method is used for evaluation of energy supply required for coplanar maneuvers. Orbit determination for two given positions of the s pacecraft The evaluation of the two burn transfer bet ween orbits 1 and 2 requires an alysis of the set of transfer orbits with further s election of the orbit providing the minimum of supplied characteristic velocity. Each transfer orbit is determined by the positions in orbits 1 and 2 and the respective times for these positions. In the scope of the non-perturbed motion the problem of orbit determination for t wo positions and the time of the trans fer bet ween them is a classical Lambert p roblem. Theorem (Euler-Lambert) For the two vectors r1, r2 of positions of a mass point in the central gravitational field and the time interval

t for transfer from position r1 to position r2 , the following equation is valid:

2


µ
where

ta

3 2

=

sin

(

sin

)

+2 n

(1)

a µ n

­ ­ ­

semi-major axis of the elliptical transfer orbit, gravitational parameter, the number of revolutions of the transfer orbit

sin

2

2

=

r1 + r2 + r2 r1 4a 0 2 ,

, sin 2

2

2 2

=

r1 + r2 4a 2 .

r2 r1

,
(2)

n n =2 + 2 n, 0 < 2 is the arc passed by the mass point, the sign "­" cor1= does not exceed , and the sign "+" to the case wh en exresponds to the case when n n ceeds . 1 and revolution number. n 2

are true anomalies corresponding to vectors r1, r2 taking into account the

Euler ­ Lambert theorem sho ws that the time for the orbital transfer of the mass point from one position to another depends only on the sum of the radii of these positions r1 + r2 , on the value of the chord connecting them r2

r1 , and on the semi-major axis of the orbit a . Equation

(1) is also called Euler-Lambert's equation. For historical referen ce1 we note that this theorem for t he parabolic orbit was first proved and published by Euler in 17435. Lamb ert's theorem, gen eralizing Euler's result was found in 17616. Albouy7 noticed that the Euler- Lambert's equation mig ht have solutions not corresponding to any transfer orbits. He suggested the following formulation of Lambert's theo rem: Four functions a , r1 + r2 , r2 r1 and t are functionally dependent. Albouy gives references to different proofs of Lambert's theorem and techniques of orbit determination for t wo positions and the time of transfer: Laplace8, Adams 9, Dziobek 10, Routh11, Plummer12, Battin13. Albouy's pedantry provides a platform fo r the applied algorithm for solving Euler-Lamb ert equation. We should search for all the solutions of the Euler-Lamb ert equation. The found solutions should be checked for compliance with the sch eme of the transfer fro m the po int r1 to the point r2 for the time interval

t . It should be noted that Albouy considers the cases when the transfer from r1 to r2 takes not more than one revolution. The evaluation of the nu mber of solutions for multi-revolution transfers is given in Ref. 14 ..
Consider the algorithm for determining the transfer trajectory for t wo given positions and the time of the transfer under the condition that the transfer is performed by not more than a given number of revolutions. This algorithm is based on the methods described in the works1,2,3 and a note of Albouy7. All the solutions of Eq. (1) are searched for. Then the found solutions are che-

3


cked for co mpliance with the scheme of the transfer from point r1 to the point r2 for time interval t . The design of the algorithm uses the following facts, 1 ,2,3: ­ we al ways have: sin

2

> 0, cos

2

>0 ;
n

the sign of sin

2

is determined by the sign of cos

2

, i . e.
n n

sign sin

2

=

1, 0 1,

< <2

(3)

cos < 0 , when the second focus of the transfer ellipse is within the elliptical sector corre2 sponding to the transfer trajectory .
Let us introduce the following notation:

sin

0

2

=

r1 + r2 + r2 4a 0
0

r1

, sin ,0

0

2
0

= 2

r1 + r2 4a .

r2

r1

,
(4)

2

2

2

Figure 1 . Segment of the Secto r Does Not Include Any Foci.

Figure 2 . Segment of the Secto r Includes the Second Focus <

Figure 3 . Segment of the Secto r Includes the Second Focus >

Figure 4 . Segment of the Secto r Includes Only the First Fo cus

4


Figure 5 . Segment of the Secto r Includes Both Foci

No w consider the following possible positions of the elliptical sector corresponding to the transfer trajectory: ­ ­ ­ ­ In t there are elliptical elliptical elliptical he case A no foci within the elliptical sector, case A, Fig ure 1, sector covers only the second focus , case B, Figures 2, 3 , sector covers only the first focus, case C, Figure 4, sector covers both foci, case D, Figure 5 . = 0 , = 0 . In the case =2 = 0 or = 0,

0

. In the case C

=

0

,

=

0

or

=

0

. In the case D

=2

0

,

=

0

.

Thus, instead of Eq . (1) we can consider four equations:

t= 1

1

µ
a
3 2

a

3 2

(

0

sin

0

(
0

0

sin

0

)

+2 n

) )

(5)

t=

µ
t= 1

(

2
3 2

(

0

sin

)(
0

0

sin

0

)

+2 n

(6)

1

µ
a
3 2

a

(

0

sin

0

+

(

sin

0

)

+2 n

) )

(7)

t=

µ

(

2

(

0

sin

0

)+(

0

sin

0

)

+2 n

(8)

Among the solutions of these equations we may find some that do not correspond to transfer trajectories, i .e. those that do not start and finish in the given points. Thus we will call these four equations an extension of the Euler-Lambert equation. The solutions of the extended EulerLambert equation should be check ed for co mpliance wi th the transfer trajectories. The selection of the numerical method for solving the equations should account for the character of the functional dependence of the right parts of the equations (5)-(8) on a . Th e functions of the right parts of the equations (5) and (7), monotonically decreas e with the increase of a for n = 0 and monotonically increase for n = 1, 2, ... . The functions of the right parts of the equations (6) and (8), increase for all values of n = 0, 1, 2, ... . Thus for solving equations (5)-(8) we can use the bisection method and the golden section method.

5


The algorithm described further was used for the task of recovery of the orbital injection scheme for the case of injection with the change of incl ination14. Algorithm fo r La mbert pro blem Input information: ­ ­ ­ ­ initial position of the spacecraft; final position of the spacecraft; time interval for the transfer; maximu m nu mber of revolutions for the transfer.

rsrc
rtrg

t N max

Output information:

N

ELM

­
k

{

k

, ik ,

, ek , pk , t

,k

}

,­

k = 1, ..., N
Wh ere
k

the number of acquired solutions (the number of reco rd s in the output array); Output array, each record contain the p arameters of the transfer orbit:

ELM

­ ­ ­ ­ ­ ­

longitude of ascending node; inclination; pericenter argu ment; eccentricity; parameter of elliptical orbit; time of passing the peri center.

ik
k

ek p
k

t

P, k

Des cription of the algorithm 1. If the vectors rsrc and rtrg are collinear, the algorithm is finished with a negative return code. Otherwise we calculate the vector m 0 =

(

0 0 mx , m0 , mz y

)

T

=

rsrc rsrc

rtrg rtrg

.

0 0 2. Calculate 2 vectors: c01 = m 0 sign mz and c02 = m 0 sign mz , 0 where sign mz

()

()

()

is the sign of the z-co mponent of vector m 0 .

The index k of the resulting array is set to zero . Items 3-5 are performed for each vector c
01

and c02 . In the description of the algorithm

these vectors will be denoted as c0 i , i = 1, 2 . Th e values calculated on their basis will have index i .as well.

6


3. Then calculate the inclination: ii = arccos c 4. The longitude of the ascending node
i
0i x

( z)
0i

0 , wh ere cz i - is the z-component of vector c

0i

is determined using conditions:

sin

i

=

c

sin ii

, cos

i

=
i

c

0i y

sin ii

(9)

5. Calculate the differen ce of true anomalies

of two given positions using the values

sin

i

and cos

i

, calculated by formulas:

cos

i

=

(

rsrc , rtrg rsrc rtrg

)

, sin

i

= c0 i , m

(

0

)

rsrc

rtrg

rsrc rtrg

.

(10)

The cycle for index i is completed. The yield is the t wo values of the difference of the true anomalies: 1 and 2. 6. Then calculat e the length s of the interval connecting the initial rsrc and final rtrg positions: . s = rsrc

rtrg .

(11)

7. Using the algorithm for solving the extended Lamb ert's equation we find two arrays of solutions:

{a {a
ond one - to
2

1, k

, ,

1, k 2, k

, ,

1, k 2, k

2, k

}, },

k = 1, ..., N

1, LEQ

k = 1, ..., N

(12)
2, LEQ

The first array corresponds to the value of the differen ce of true anomalies .

1

, and the s ec-

8. Then we calculate the unit vector corresponding to vector rsrc
0 rsrc =

(

0 xsrc ,

0 ysrc , z

0 src

)

T

=

rsrc rsrc

(13)

9. The cycle for j = 1, ..., N

1, LEQ

+N

2, LEQ

performs op erations of items 10-16 .

10. Five values are generated:

ij =

i1 , if j

N

1, LEQ j 1, LEQ

i2 , if j > N

=

1 2

, if j

N

1, LEQ 1, LEQ

, if j > N

aj =

a1,

j jN
1, LEQ

, if j
+1

N

1, LEQ

a2,

, if j > N

1, LEQ

7


j

=

1, j 2, j N
1, LEQ

, if j
+1

N

1, LEQ j

, if j > N

=

1, j 2, j N
1, LEQ

, if j
+1

N

1, LEQ

1, LEQ

, if j > N

1, LEQ

11. Then calculate the argument of latitude:

u=

0 arccos( xsrc cos

j

0 + ysrc sin j

j

),
j

if z

0 src 0 src

0 <0
and

2

0 arccos( xsrc cos

0 + ysrc sin

), if z

(14) with the values

12. Then we calculate the parameters connecting the auxiliary angles of eccentric anomaly in the initial and final points.

E2
e cos E2 p1 = ecE

E1 = E2
=1

m1

=

j

j

2
1 cos E2
, E2 p1 =
m1

(15)

rsrc + rtrg 2a
j

2 p1

(16)
m1

e sin E2 p1 = esE

2 p1

=

rsrc 2a
j

rtrg

1 sin E2

E1 + E2 2

(17)

Using the values of the sine and cosine of the h alf-sum of eccentric ano malies we can calculate the v alue of this half-sum E2 p1 . For this purpose we can use the function atan2 ( ). The values of the eccentric anomalies for the initial E1 and the final E2 points are calculated using formulas:

E1 = E2

p1

E2

m1

(18) (19)

E2 = E2 p1 + E2

m1

13. Then we calculat e: the value of the eccentricity

ej = cos

1 E1 2

rsrc + rtrg 2a E2 E + E2 cos 1 2
(20)

14. the true ano maly fo r the initial point:

= 2arctg

1+ e 1e

j j

tg

E1 2

(21)

8


15. the argument of the pericenter:
j

=u

(22)

16. time of p assing the pericenter:

t

P, j

=

a

3 j

µ

(

E1 e sin E1

)

(23)

17. parameter o f the elliptical orbit:

pj = aj 1 e

(

2 j

)
r
j ,1

(24)

18. Using the determined elements we calculate the state vectors for the times t1 = 0 and

t2 = t :

r

j ,1

and r

j ,2

. Then we calculate the residual: d j = rsrc

+ rtrg

r

j ,2

. If d j is

smaller than the set threshold value the solution is saved, otherwise - discarded. In case the solution is saved we mak e the following assignments:

k = k + 1,

k

=

j

, ik = i j ,

k

=

j

, ek = e j , pk = p j , t

P ,k

=t

P, j

19. After the completion of the cycle for j index k is used for forming the length of the output array: N ELM = k . Algorithm fo r the Extended La mbert Equation The algorithm finds the solutions of the extended Lambert's equation (1). The algorithm generates the array containing the values of the semi-majo r axis and auxiliary angles and . Th e true solutions of Lambert's equ ation, i.e. the solutions corresponding to certain transfer o rbit are among the elements of this array. The s earch for the t rue solutions requires an additional check which is p erformed by the algorithm of higher lev el. Input information:

rsrc

­ ­ ­ ­ ­

Magnitude of the vector defining the initial position of the spacecraft , Magnitude of the vector defining the final position of t he spacecraft, time interval o f the transfer, maximu m of the nu mber of revolutions for the transfer, difference of true anomalies.

rtrg t N max

Output data:

N

{

LEQ k

­

the number of elements in the array, determined values.

ak ,

,

k

,N

k

}

k = 1, ..., N

LEQ

­

We set the initial value of index k = 0 .

9


The operations 1-3 are performed in the cycle for N from the minimum possible value of the semi-major axis a = and the angle 0 <
0

0to N max 1 . We calculate
of the transfer trajectory

(

)

rsrc + rtrg + s 4

<

, determined by the equality:

sin

0

2

=

rsrc + rtrg 4a

s

(25)

The minimu m value o f the semi-major axis can be attained only for the case boundary elliptical sector. In this case the segment connecting the ends of radius- vectors r1 and r2 , includes the second fo cus of the ellipse, the angle = . If , we check the equ ality:

µ t=a

3 2

(

0

sin

0

)

+2 N
=,=
k 0

(26)

If the equality is satisfied we save in the output array the values: a, ing this we add 1 to the index k and make the assignments: ak = a, After that we return to the beginning of the cycle with new value o f N . If

, N . For do0

=,

k

=

, Nk = N .

<

2 , we check the equality:

µ t=a

3 2

+

(

0

sin

0

)

+2 N
k

(27)

If the equality is satisfied we make the assignments: k = k + 1 , ak = a ,

=

,

k

=

0

,

N = N k and we go to the beginning of the cycle with new v alue of N .
If none of the above mentioned equalities are satisfied we go to the next item 2 . 1.

+ a to amax we search for the zeroes of 4 two functions whose shape depends on the value of th e angle . Here the parameter a
is of the ord er of unit meters, and amax does not exceed 300 ,000 km. If , we search fo r the zero es of functions:
1

Within the interval fro m amin =

rsrc + rtrg + s

( a)

= µt a
3 2

3 2

0

sin

0

(
0

0

sin 0 + 2 N ,
0

)

(28)

2

( a)

= µt a

2

(

0

sin

)(

sin

0

)

+2 N .

If

<

2 , we search for the zeroes of functions:

10


3

()

a= µt a
3 2

3 2

0

sin

0

+

(

0

sin

0

)

+2 N ,
(29)
0

4

( a)
sin

= µt a

2

(

0

sin

0

)+(
=

0

sin

)

+2 N .

Wh ere
0

2

=

rsrc + rtrg + s 4a

, sin

0

rsrc + rtrg 4a

s

2

. we find the value a

(30) ,

Thus, using the bisection algorithm within the interval amin , amax for which the function
k

=0

is equal to zero .
0

3. Then we calculate the angles 0 <

< , 0<
, sin
0

0

<

, determined by the equ alities:

sin

0

2

=

rsrc + rtrg + s 4a
=0

2

=

rsrc + rtrg 4a
=0

s

. and :

(31)

Dep ending on the nu mber of the function k we calculate the angles

k
1 2 3 4 ing values of
k

=0

=
0

0

2 =
2
0

= = =

0 0 0 =0

0

In the output array we sav e the found value of the semi-major axis a

and the correspond=0

, , N . For this we make the assignments : k = k + 1, ak = a

,

k

=,

= , Nk = N .

4. After completion of the cycle the index k is used for forming the length of the output array: N LEQ = k . Evaluation o f the energy costs for co planar tra nsfers. Geo metrical metho d For maneuver characterization problems we al ways k now the state v ectors for the spacecraft before and after the maneuver. In the scope o f spacecraft control theory these p roblems are called rendezvous problems. Thus we should use the techniques for solving rendezvous problem as the basis for the maneuver characterization problem. In the case of the non- p erturbed motion this is the Lambert problem. Ho wever, for ev aluation of the results we should compare the obtained energy costs with the case for which we have the parameters of the targ et orbit, but the position in this orbit is not defined. Wh en we control the spacecraft for reaching certain position for a given time we perfo rm phasing maneuvers which provides such conditions for the orbital transfer that the energy cost for the orbital transfer with the set position in the target orbit does not signifi-

11


cantly differ from the orbital transfer with an arbitrary position of the spacecraft in the target orbit after the transfer man euver. Let us consider the geo metrical method of the transfer fro m one coplanar o rbit to another. Consider the family of ellipses with the focus in the origin of coordinates. This family of ellipses is described by the parameters l and f which are determin ed by the following relationships

c=

f l c , a= , e= = , 2 aa 2

f

(32)

where c ­ half distance bet ween the foci a ­ semi-major axis e ­ eccentricity If the p ericenter is placed to the left of the coordinate o rigin then f > 0 , otherwise f < 0 The distance to the spacecraft and the velocity at apocenter r , V have the following relationships with l and f :

.

and pericenter r , V

r= r=

l+ f 2 l 2 f

,V= ,V=

2µ l f , l l+ f 2µ l + f , llf

(33)

where µ is the gravitational constant. Each orbit is represented by a point in the semi-plane l , f , l > 0 . If f > 0 , we have apocenter right from the origin of coordinates, if f < 0 - pericenter. We will consider applying the burn pulses only at pericenter and apocenter. Let us consider the orbital transfer from the point l1 , f1 to the point l2 , f 2 plied at the apsidal point right-side from the origin of coordinates, the orbital the distance to this apsidal point. For f1 > 0, f 2 > 0 - this is an orbital tran apocenter will remain right-side from the origin of coordinates before and after . If tran sfer the the burn is apsfer must keep for which the burn , thus (34)

l1 + f 2

1

=

l2 + f 2

2

The case f1 > 0, f 2 < 0 correspond to the orbital transfer for which we will have the pericenter right-side fro m the origin of coordinates after the bu rn is applied. Thus:

l1 + f1 2

=

l2 2

f

2

=

l2 + f 2

2

(35)

Consideration of the cases f1 0, f 2 < 0 and f1 0, f 2 > 0 , yields that the application of the burn on the right-side of the o rigin of coordinates keeps the value l + f .

12


Similar considerations on the cases when the bu rn is applied left-side from the o rigin of coordinates will result in the conclusion that for these cases the value l + f is kept. Thus the broken line in the semi-plan e l , f , l > 0 , corresponding to the sequence of orbital transfers will consist of segments with inclination angle tangents -1 and +1 (Fig. 6). The value -1 for the inclination angle tangent corresponds to the burn applied right-side from the origin of coordinates and the value +1 ­ corresponds to the left-side application of b urn. The v elocities in the apsidal points (in the apocenter or pericenter) right-side Vr and left-side Vl from the origin of coordinates are calculated using formulas:

Vr = Vr = 2 µ ( l2 f2 )

2µ(l f ) , Vl = l (l + f ) 2 µ ( l1 f1 ) , Vl =

2µ(l + f ) . l (l f ) 2 µ ( l2 + f 2 ) l2 ( l2 f2 ) 2 µ ( l1 + f1 ) l2 ( l1 f1 ) .

(36)

l2 ( l2 + f 2 )

l1 ( l1 + f1 )

(37)

For the orbital transfer from the point l1 , f1 to the point l2 , f 2 the values of the transfer impulses Vr (when the burn is applied right-side fro m the origin of coordinates) and Vl (when the burn is applied left-side fro m the origin of coordinates) are calculated using the formulas:

Figure 6 . Orbital Transfers D epict ed In the SemiPlane l , f , l > 0 . Fo r the Transf er From the Point l1 , f1 To the Po int l2 , f 2 The Burn is Applied R ight-Side From the Origin Of Coordinates, For the Transfer Fro m l2 , f 2 To the Point

Figure 7 . Hohmann's Transf er Plotted In the Semi-Plane l , f , l > 0 . The Upper Broken Line Corresponds To the First Burn Applied LeftSide From the Orig in of Coordinat es, the Low er Line ­ To the Right-Side Application.

l3 , f3 ­ To the Left Side
Let us consider the Hohmann's transfer fro m the circular orbit with radius l1 to the circular orbit with radius l2 (Figure 7). The elliptical transfer orbit will have the p arameters:

l3 =

l1 + l2 2

, f3 = ±

l2 2

l1

. The sign « +» correspond to the burn applied to the left-side of the

13


origin of coordinates and the sign «-» - to the right-side application. Using (37) for l2 will obtain the consumption of the characteristic velocit y for the Hohmann's transfer:

> l1 , we

VH =

2µ C l1

()

,

C

()

=

2 1+

1+

1

(

2

+1

)

,

(38)

where

=

l2 l1

.

The graph o f the function C

()

is shown in Figure 8.

Let us consider the three burn bi-elliptical transfer from the circular orbit with radius l1 to the circular o rbit with radius l2 (Figure 9). The first burn is applied left-side fro m the origin of coordinates and makes the semi-major axis of the first transfer orbit where times greater, i.e. l

p1

= l, 1

> 1 . Parameters of intermediate orbits l p1 , f p1 , l p 2 , f

p2

are connected by the followin g

relationships:

l1 2

=

l

p1

f 2

p1

,

l

p1

+f 2

p1

=

l

p2

+f 2

p2

,

l

p2

f 2

p2

=

l2 2

.

(39)

Figure 8 . Graph Of the Funct ion C

()

.

Figure 9 . The Three Burn Bi- Elliptical Transfer From the Circula r Orbit With Radius l1 To the Circula r Orbit With Radius l2 and the Hohmann's Transf er. Hohmann' s Transf er Corresponds To the Broken Line Including the Points: (l1,0) - > ( l3,f3) -> (l2,0). The Three Burn Bi- Elliptical Transfer Co rresponds To the Broken Line Including the Po ints: (l1,0) -> (lp1,fp1) - > (lp2,fp2) -> ( l2,0).

Solving the system of equations (39) with respect to f p1 , l

p2

,f

p2

, we will have:

14


f

p1

=

(

1 l1, l

)

p2

=

(

2

1 l1 + l2 2

)

,f

p2

=

(

2

1 l1 2

)

l2

(40)

The burn impulses for the first, the second and the third segments of the bi-elliptical transfer are calculated using formulas:

V1 =

2µ l1

2

1

1,

V2 =

2µ 1 l1 2 1

(

2l2 2

1 l1 + l

)

1
2

,

(41)

V3 =

2µ 1 l2

(

22 2

(

1 l1 + l2

)

1 l1

)

The total v for the three burn bi-elliptical transfer is calculated by the formula:

V3 B =
For

V1 +

V2 +

V3

(42)

= l1 + l2 /2 the bi-elliptical three burn transfer degenerates into the Hohmann's trans-

(

)

fer and the v alue of

V3 B beco mes equal to

VH .

The non-dimensional values l1 , l2 - mean the ratios of the radii of the o rbits to mean Earth radius. For transition to the dimensional value of the magnitude of the ch aract eristic velocity we should multiply the non-dimensional value by the value of the first escap e velocity

µ / RE = 7.905365716 km/s .
We will co mpare the ch aracteristic velocity for the th ree burn bi-elliptical transfer with that for l1 = 1, l2 = 2 for the Hoh mann's transfer. Figure 10 sho ws the dependence of V3 B from with the characteristic velocity for the Hohmann's tran sfer at the b ackground which are presented by the straight line parallel to the abscissa axis. For = 1.5 the three burn transfer degen erates into two burn Hoh mann's transfer with energy consumption of 0.67884128 in non-dimensional greater or lower than 1 .5 the consumption for the three burn transfer values. Fo r the values of is greater than that fo r the Hohmann . Let us now set

= 1 + l2 and compare the characteristic velocity for the three burn and

Hoh mann's transfer from the circular orbit with radius l1 = 1 to the orbit with radius l2 . Figure 11 shows the differen ce of the ch aract eristic velocity for the three burn and Hohmann 's transfer. Within the interval of l2 values from 13 to 14 the difference in the characteristic velocity changes from 0 .0029183931 to -0.0006411714, i.e., it changes sign. Thus for the value o f l2 = 14 we have less change in characteristic velocity for the three burn bi-elliptical transfer than for the Hoh mann transfer.

15


Figure 10. Co mparing The Chara cterist ic Velocity For the Three Burn Bi-Ellipt ical and Hohmann's Transf er For the Radii of 1 And 3 For Init ial and The Target Orbit Respect ively. The Abscissa Ax is Present s - The Ratio Of the Pericenter D istance of the First Transfer Orbit To the Radius of t he Init ial Orbit . The Ordinate Axis Presents the C haract eristic Velocity Consumption In Non-Dimensiona l Characterist ic Velocity Consumpt ion For the Hohmann's Transfer.

Figure 11. The Diff erence of C haract eristic Velocity (In Non-D imensional Values) . Fo r the Three Burn Bi- Elliptical and Hohmann's Transfer As Funct ion of l2 for l1 = 1 , = 1 + l2

The gen eral result is known [4]. For the ratio of orbital radiuses l2 / l1 < 11.94 the characteristic velocity for the Hohmann transfer is smaller than for the three burn bi-elliptical transfer. For the ratio of the radii 11.94 < l / l < 15.58 there exists the value of , for which the consump21 tion for the three burn transfer is smaller than for the Hoh mann's one. For the ratio of the orbital radii l / l > 15.58 the characteristic velocity for the three burn transfer is smaller than for the 21 Hoh mann transfer for any value of > 1 . Ho wever, the increase in the characteristic velocity by using the three burn bi-elliptical transfer instead of Ho hmann's one does not exceed 8 %, and the time required for the transfer increases sev eral ti mes. The above considerations lead to the following conclusions. The evaluation of the characteristic velocity for the Hohmann 's transfer is the basic val ue for checking the reliability of other possible maneuvering pro files with regard to energy cost. The geo metrical technique for analysis of orbital transfers is an efficient tool for analysis of maneuvering profiles without changes of the orbital plane. The gen eral result is known [4]. For the ratio of orbital radiuses l / l < 11.94 21 tic velocity for the Hohmann transfer is smaller than for the three burn bi-elliptic the ratio of the radiuses 11.94 < l / l < 15.58 there exists the value of , for 21 sumption for the three burn transfer is smaller than fo r the Hoh mann's one. For the characterisal transfer. For which the conthe ratio of the

16


orbital radii l / l > 15.58 the characteristic velocity for the three burn transfer is smaller than 21 for the Hohmann transfer for any value of > 1 . Ho wever, the increase in the characteristic velocity by using the three burn bi-elliptical transfer instead of Hoh mann's one does not exceed 8%, and the time required fo r the transfer increases sev eral times.

Figure 12. D ifference Of Cha racteristic V elocity Consumption (In Non-Dimensiona l Va lues) Fo r the Three Burn Bi- Elliptical and Hohmann's Transf er As Function of l1 and l2 for l1 = 1 , = 1 + l2

Figure 12. presents the difference in the characteri stic velocity for the three burn and = l1 + l2 . Hoh mann's transfers for 1 l1 1.1, 5 l2 20, The above considerations lead to the following conclusions. The evaluation of the characteristic velocity for the Hohmann 's transfer is the basic val ue for checking the reliability of other possible maneuvering pro files with regard to energy cost. The geo metrical technique for analysis of orbital transfers is an efficient tool for analysis of maneuvering profiles without changes of the orbital plane. ALGORITHM FOR EVALUATI ON OF ONE BURN MANEUVER The ch ange of orbital parameters caused by thrust forces, in the case the orbital parameters (TLE) before and aft er the burn are available can be i nterpreted as a on e burn maneuver. This is

17


the most frequent case. If the determination of the orbital parameters occurred more frequently than the maneuvers , all the ch anges in the orbits could be interpreted as one burn man euvers. If a one burn maneuver occurred the trajectories before and after the burn must intersect in the point where the thrust was applied. Thus the an alysis of possible causes of the ch ange of the o rbital parameters we should search for the point corres ponding to the minimum distance bet ween the trajectories before and after the thrust. The resulting vector of residuals between orbital positions before and after the burn should be compared wit h the error of position prediction based on the orbital parameters before and after the burn. Algorithm fo r the search of the mi ni mum dista nce between the orbits Wh en we search for the minimu m distance bet ween t he orbits we know the time interval for the search. This interv al coincides with the time interval within which the burn was applied. It should be noted that in case we have available the TLE before the thrust and the TLE after the thrust the epochs of these TLEs should not be identified with the interval of maneuv ering. The TLE epoch normally refers to the time of nodal crossing and not to the times of the last measurements of the observation interval for whi ch the TLE h ave been generated. Thus in case we hav e the TLE, the interval for the search o f the minimum dis tance should be extended by on e period to the left from the epoch of TLE before the burn and by one period to the right from the epoch of TLE after the burn . Let us introduce the notation:

t1, t

2

­

r1 t =
v
1

()

r2
v

2

( t () = ( t () = ( t () = (

x1 t ,

()

y1 t , z1 t
y1

()

()
z1

)

T

­

Vx1 t , V

()

()

t, V

()
t

x2 t ,

()

y2 t , z 2 t
y2

()

()
z2

)

)

Boundaries fo r the search of the minimum distance bet ween the orbits; Position vector of the spacecraft based on orbital data before the burn; Velocity v ector of the spacecraft b ased on orbital data before the burn; Position vector of the spacecraft based on orbital data after the burn; Velocity v ector of the spacecraft b ased on orbital data after the burn.

T

­ ­

T

Vx 2 t , V

()

t ()

,V

t ()

)

T

­

The task is to find the time t , for which the v alue r2 t

()

r1 t

()

reaches its minimum.

r2 t

()

For the search of the minimum distance we will use the fact that the derivative of the function

r1 t

()

must ch ange sign. In the cas e when for the entire interval t1, t
2

2

the derivative o f

the considered function does not ch ange sign, the minimum distance corresponds to one of the boundary values of the interv al t1, t by the function: . The d erivative of the function r2 t

()

r1 t
V

()
)

is calculated

d r dt 2

r1 =

(

x

2

x1 V

)(

x2

V

x1
2

)+(

y

2

y1 V
2

(

x

x1 + y

)(
2

)(

y2
2

V

y1 1

y1 + z

)(

)+(
2

z

2

z1 V

)(

z2

z1

z1

)

2

(43)

18


The s earch for the points where the d erivative of the function r2 t

()

r1 t

()

changes sign is

performed by scanning of the interval of the search. Howev er, the scanning should be performed not by the time as variable, but using the true anomaly with permanent step h , equal, for example, to 1° or 10°. Consider the transition from the time of the previous step t p to the time of the current step of scanning t . For the time t p using orbital data before the maneuver we calculate the osculating orbital elements:

M

p

­ ­ ­ ­

mean ano maly; true ano maly; eccentricity; mean motion.

p

e n

p p

The transition to the time t of the current step of s canning is performed using formulas:

E = 2arctg

1e 1+ e

p p

tg

p

+h 2

M=

E e p sin E ,

if E e p sin E

M p,

E e p sin E + 2 , if E e p sin E < M p ,

(44)

t = tp +
The scanning is finished wh en t > t2 .

M n

M
p

p

For each step of s canning we calculate the d erivative of the distance b et ween the orbits:

Dt=
using Eq. (43). If the condition:

()

d rt dt 2

t () r ()
1

(44)

Dt

t ( ) D ()
p

0

(45)

is satisfied the derivative ch anges its sign within the interval t p , t . In this case for the time

tm = tp + t /2 we calculate the value of the derivative of the distance between orbits D t

(

)

( m)

.

19


If the values D p = D t zero, using the formula:

()
p

, Dm = D t

( m)

, D = D t are monotonic, i .e. D p < Dm < D

()

or

D p > Dm > D we calculate the time t0 , when the derivative of the distance between orbits is

t0 = t

p

(

Dm D D
p

Dm

)(

D

p

D

)

+t

m

(

D Dm D
p

p

)(

D Dm D

)(

+t

Dm D D Dm

)(

DD

p

)

.

(46)

If the monotonic condition for the values o f D p , Dm , D is not satisfied, we set: t0 = tm . Calculation of the i mpulse for the scheme of the one burn maneuver Further we perform the calculation of the consumption of characteristic velocity for the maneuver and the vectors of relative positions and velocities in the RNB coordinate frame. For this purpose for the time t0 we calculate the state vectors

(

x2 ,

y2 , z2 , Vx2 , Vy2 , Vz

2

)

(

x1 ,

y1 , z1 , Vx1 , Vy1 , Vz1

)

T

,

T

using initial conditions before and after the maneuver.

Then we calculate the consumption of the characteristic velocity for the maneuver using the formula:

DV =

(V

x2

Vx1 + Vy2 Vy1 + Vz2 Vz1

)(
2

)(
2

)

2

(47)
RNB

On the basis of the state vector befo re the maneuver we calculate the matrix M transformation to the RNB coordinate frame: M
RNB

of the

=

er , e n , e

T b

, where

er = eb = r1 r1
r

r1 r1 V1 V1

­ ­ ­

unit vector directed to the point r`1 =

(

x1 ,

y1 , z1

)

T

,

unit vector orthogonal to the orbital plan e, V1 =

(

Vx1 , V y1 , Vz

1

)

T

,

en = e

eb .

The v ectors of the relative positions and velocities for the RNB coordinate frame are calculated by the formulas:

x DR
RNB

2 2 2

x1 y1 z1
(48)

=M

RNB

y z

Vx DVRNB = M
RNB

2

Vx V Vz

1

V

y2 2

y1 1

(49)

Vz

20


Reliability criteria for the obtained evaluation of the impulse The major criterion to be used for assessment of the reliability of the obtained evaluation of the impulse is the distance bet ween the positions of the spacecraft in the first and the s econd orbit for the time of minimum distance. This distance is calculated by the fo rmula:

DR =

(

x

2

x1 + y

)(
2

2

y1 + z

)(
2

2

z1

)

2

(50)

If the value DR exceeds the threshold determined by the accuracy of the orbital data we can consider that by the time o f minimum distance the orbits do not intersect and the obtained estimate is not reliable. Wh en DR is smaller than the threshold we perform additional analysis including comparing the value of DV with the characteristic velocity required for the Hohmann transfer or inclination changing maneuver. If the orbital planes before and after the man euver coincide, a coplanar maneuver has occurred . En ergy cost of this maneuver should be comparable with the cost of one of the impulses of the Hohmann transfer. If the orbital plane has changed the value of DV should be co mpared with the characteristic velocity required for the inclination changing maneuver. We should note one very important empirical criterion for the impulses with modules not exceeding 1 m/s. If the maximu m (in absolute value) co mponent of the vector DR RNB is the co mponent corresponding to the direction of vector N , and the maximu m (in absolute value) component of the vector DVRNB - is the component corresponding to the direction of vector R , the obtained estimate of the maneuver is not reliable and the difference bet ween the orbital parameters is caused by the errors of orbit determination. Algorithm fo r the calculation of the length of the i mpulse The length of the impulse can be calculated knowing the characteristic velocity, if the mass of the spacecraft b efore the burn, the thrust and the specific impulse of the engine are kno wn. Under the conditions when these parameters are not availabl e the length of the impulse can be determined on the basis of the following empirical table wh ich presents the acceleration as function of the characteristic velocity. The table is based on the simple consideration that small impulses would not be generated by high-powered engines and the big ones ­ by the lo w thrust engines. Characteristic velocity, m/s 0-5 5-70 70-150 150-1000 1000-... Dep ending on DV - the estimate of the characte Acceleration, m/s2 0.1 0.2 0.5 2 10 ristic velocity of the maneuver, the mean ac-

celeration a is determined. The length of the impulse time of the start of the man euver t formulas:
beg_imp

t

imp

is calculated as

and the time t

end_imp

DV . The a of its end are calculated using the t
imp

=

t

beg_imp

=t

t
0

imp

2

,t

end_imp

= t0 +

t

imp

2

(51)

21


where t0 is the time corresponding to the minimum distance bet ween o rbits (the time of the burn in the scope o f impulse understanding) ALGORITHM FOR EVALUATI ON OF TWO BURN MANEUVER , t E 2 , t E1 < t E 2 and the i nterval of possible maneuvering tm1 , tm2 are known. The orbit corresponding to the time t E1 , will be called the initial one and the orbit corresponding to the time t E 2 - the target one. The position vector of the spacecraft as a function of time will be called the trajectory. Corresponding to the initial and target orbit we will differ the initial and the target trajectories. If within the interval tm1 , tm2 we cannot find the time of closest approach o f the initial and target trajectories we will consider that between the times tm1 and tm2 more than one burn h as been performed . Then we will consider the algorithm for the evaluation of the impulses and the times of their application under the assumption that a two burn maneuver has been performed resu lting in the transfer from the initial orbit to the target one. Assume the orbital parameters for times t
E1

First we will consider the variant when the orbital parameters are represented by position r1 , r2 and velocity v1 , v 2 vectors for the times t E1 , t E 2 . The task of ev aluating the impulses and the times of their application is equivalent to the task of determining the transfer orbit between certain positions in the initial and target trajectories respectively. Selecting the times t , t t ,t , t < t we can obtain the infi nite set of such 12 m1 m2 1 2 orbits. The desired orbit should be selected to satisfy the condition of the minimum consumption of characteristic velocity for transfer. Selection of the time t1 defines the position of the spacecraft in the initial trajectory and the selection of the time t2 defines the position of the spacecraft in the target trajectory. Thus for the search of the transfer orbit we have t wo positions and the respective times. So the search of the transfer orbit can be considered the solving of the boundary problem for the differential equation describing the motion of the spacecraft. Since here we do not deal with the task of spacecraft flight control and consider the task of evaluation of the maneuvers which have been performed already, we can look for approximate estimates. Ho wever, we should mention that the obtained approximate estimate is the initial approximation fo r ob taining the precise solution. The considered problem is solved according to the following scheme. We scan the set of transfer orbits and look fo r the orbit corresponding to the minimum change in characteristic velocity. The transfer o rbits are generated approxi mately using two techniques. The technique fo r generating the transfer orbits is selected depending on the condition on the value of the impulses required for transition to the transfer orbit from the initial one and for transition from the transfer orbit to the target one. We consider the following cases: ­ the values of the impulses are smaller than the given threshold and we can use linearity in the vicinity of the reference trajectory; ­ the occurred ch ange of the orbital parameters corresp ond to the change in characteristic velocity exceeding the set threshold.

(

)

22


Wh en we h ave the opportunity to use the linearity in the vicinity of the reference trajectory then the basic tool for the approximate gen eration of the transfer orbit is the matrix of partial derivatives of the components of the state vector for the time t2 with respect to the components of the state vector for the time t1 , calculated along the initial orbit. If we can not use the linearity then the basic tool for generating the transfer trajectories is Lambert's problem, which can find the transfer orbit between two positions for the given time interval of the transfer. The interval of the possible application of the first o r the second burn can be limited for the following cas es: the initial orbit is an eccentric one and the target orbit is near-circular, the inclinations of the orbits coincide; the initial orbit is a near-circular one and the target o rbit is eccentric, the inclinations coincide; the inclinations of the initial and target orbits are different and the semi-major axes are virtually the same. Let us consider the major idea of the algorithm based on linearity. Let for the time t1 we kno w the state vector of initial orbit rS t1 , v
T ( ( ) S (t1 ))

and the matrix of partial derivatives of the com-

ponents of initial orbit for the time t2 with respect to the components of the state vector for the time t1 :

(

t2 , t1 =

)

(rS (t2 ) , v S (t2 ) ) (rS (t1 ) , v S (t1 ) )
V1 , and at the time t2 - the impulse

(52)

If at the time t1 we apply the impulse

V2 , the state

vector for the time t2 can be approximately represented i n the shape:

rS t v
S

( 2) ( t2 )

+

(

t2 , t1

)

0 V1

+

0 V2

(53)

If at the time t2 the spacecraft transferred to the target orbit, the approximate equality is valid:

rT t v
T

( 2) ( t2 )

rS t v
S

( 2) ( t2 )

+

(

t2 , t1

)

0 V1

+

0 V2

(54)

We will consider this equality as an equation for the determination of the vectors ate them and det ermine the change in characteristic velocity:

V1 and

V2 . Thus, assuming that the thrusts have b een p erformed by the times t1 and t2 , we can evaluV= V1 + V2 . The times of
the thrusts are d etermined by scanning the interval of their possible application under the condition of the mini mum consumption of characteristic velo city.

23


Algorithm fo r evalua tion of two burn maneuv ers using linearity Input information:

(
t

t E1 , rEI , v
t1L , t
E2
2L

(

)

EI

)T )T

­ ­ ­ ­

the time and the state vector of initial orbit for this time; interval of possible thrusts for transition from initial orbit to the transfer on e; the time and the state vector of the target orbit for this time; interval of possible thrusts for transition from transfer orbit to the target one;

(
t t

, rE 2 , v
2R

(

t1R , t

)

E2

Output information:
I1

­
1

evaluation of the ti me o f the first thrust; evaluation of the first thrust; evaluation of the ti me o f the s econd thrust; evaluation of the s econd thrust. state vector o f the initial orbit for the time t ; state vector o f the t arget orbit for the time t ;

VI
I2

­ ­ ­ ­ ­

VI 2 Designations

t t (r () , v ()) t t (r () , v ())
S S T T

Algorithm 1. For the time interval t1R , t

(

2R

)

using the orbital parameters of the target orbit we calculate the

table of the state vectors with fixed step: tT ble comp rises N
T

{

,i

, rT tT ,i , vT tT

( ) ( ,i )}

. We will co nsider that this ta-

elements.
T

2. For each time tT ,i , i = 1, ..., N

using the orbital parameters of the initial orbit we calculate the
,i

matrix of partial derivatives of the co mponents of the s tate vector for the time tT the components o f the state vector for the time t

with respect to

(

tT ,i , t

E1

(r (t ) , v (t ) ) ) = (rS (tT ,i ) , v S (tT ,i ) ) S E1 S E1
rT tT v
T

E1

:

(55)

and the vector of residuals between the state vectors i n the target and initial orbits for the time tT ,i :

r tT v

( ,i ) (tT ,i )

=

( ,i ) (tT ,i )

rS tT v
S

( ,i ) (tT ,i )

(56)

24


3. Within the interval t1L , t

(

2L

)

using parameters of the initial orbit we calculate with permanent

step the state vectors and the matrix of partial derivatives of the state vector for the initial time with respect to the state vector for the current time. The results are saved in the table. We will consider that the table comprises N L elements. In this table for each time t L, j , j = 1, ..., N L we have:

rS t

( ) , v (t )
L, j S L, j

and

(

t E1 , t

L, j

r (t ) , v (t ) ) = (r t , v t ) ( ( ) ( ))
S E1 S S E1 S L, j L, j
T

(57)

4. Then we p erform the enu merative calculations for i = 1, ..., N of the enu meration we calculate the matrix:

and j = 1, ..., N L . For each step (58)

(

tT ,i , t

L, j

)= (

tT ,i , t

E1

)(
0 V1

t E1 , t

L, j

)
0 V2

and solve the equation:

r tT v
If the matrix

( ,i ) (tT ,i )

=

(

tT ,i , t

L, j

)

+

(59)

(

tT ,i , t

L, j

)

is represented in the shap e of block matri x:

(
r tT

tT ,i , t

L, j

)

=

11 21

12 22

(60)

the considered equation splits into two equations:

( ,i )
=

=
1

11

V1 , v tT ,i =

()
V2

22

V1 + V2

(61)

We can find the solution of this system of t wo equation s using formulas:

V1 t
ules of impulses:

()
L, j

11

r tT

( ,i )

,

(tT ,i )

= v tT

( ,i )

22

V1

(62)

The result of the enumerative search in the pair of indi ces im

, jm , for which the sum of the mod-

V1 + V2 reaches the minimum.
t t
I1

5. Then the result is gen erated:

= t L, j , m = tT ,i , m

VI 1 = V1 t VI 2 =

I 21

( L, j ) , V2 (tT ,i )

.

(63)

Algorithm for evaluation o f the two burn maneuv ers in the case of significa nt chang e of orbital pa ra meters In the cases when the man euvering changes the o rbital parameters to the extent that we can not use linearity, the thrusts have such values that making their evaluation we can neglect the non-central features of the Earth gravitational field. In this case we can use the Lambert's problem. If a precise estimate is required the solution obtained for Lamb ert's problem is the initial approximation for the iterative scheme. The algorithm for evaluation of t wo burn maneuvers us-

25


ing Lambert's problem is similar to the procedure considered above. The difference is that instead of solving the Eq . (54) we solve Lambert's problem. Using the evaluation of the ti me of the transfer fro m the initial to the target orbit for reducing the enumerative search 1. Minimum time for the transfer is calculated by the fo rmula:
3 atrn

ttrn =
where

µ

(64)

atrn =

a1 + a2 2

-

semi-major axis of the transfer orbit, semi-major axis of orbit 1, semi-major axis of orbit 2.

a1 a2

2. Co rrection of the interval for the search o f maneuvers. Denote the initial interval for the search as t1 , t2 . If t2 t1 < 2ttrn , then the correction of the boundaries of the interval is needed . Correction of the boundaries of the search interv al is performed using the following procedure. We find the time tmin , corresponding to the minimu m distance bet ween the orbits. It is expedient to use for t
min

the value determined by evaluation of the one burn maneuvers. If t1 > t

min

ttrn , we

change t1 ; and set it to tmin ttrn . If t2 < tmin + ttrn , we hav e to change t2 . t2 is set to tmin + ttrn . 3. For redu ction of the enumerative search we should analyze only such times for the 1st and the 2nd thrusts for which the time interval b etween them ex ceeds ttrn . Using the features of the i nitial and targ et o rbits for reduction of the enumerative search The intervals fo r the enumerative search can be reduced in special cases. These are the following cases: ­ transfer from the eccentric orbit to the n ear-circular on e without change of the inclination (circling the orbit); ­ transfer fro m the near-circular orbit to the eccentric one without change of the inclination; ­ change of the inclination without significant change of the semi-major axis (less than 1000 km). Transfer from the eccentric orbit to the n ear-circular o ne without change of the inclination The 1st thrust is looked for in the vicinity:

180
where
1

0 1

1800 +

(65)

is true anomaly of the 1st orbit, adjusted to the i nterval 0º - 360º.

Transfer from the n ear-circular orbit to the eccentric o ne The 2nd thrust is looked for in the vicinity (of the apocenter):

26


180
where
2

0 2

1800 +

(66)

is true anomaly of the orbit 2, adjusted to the interval 0º - 360º.

Change of the inclination without significant change of the semi-major axis (less than 1000 km). The 1st thrust is looked for within the following interval s:

0 u1 180
0 0

u, u u1 1800 + u , u u1 360 ,
0

(67)

360

where u1 is the argument of l atitude of the sp acecraft in the orbit 1. Parameter determines the vicinity of the apocenter and parameter u -- the vicinity of the node. The default values of these parameters: 90º. Using the interactive mode these values can be changed . The algorithm for evalua tion of two burn ma neuv ers for the case of simultaneous cha nge of semi-ma jor axis a nd inclination. Input information:

t
t

E1
E2

, rE1 , v

(

E1

)T )T

­ ­

the time and the state vector of the initial orbit for this t ime; the time and the state vector of the target orbit for this time.

, rE 2 , v

(

E2

Output information:

ta1 , Va

1 2

­ ­

the time and the components o f the 1 st thrust, the time and the components o f the 2
nd

ta 2 , Va

thrust.

Algorithm 1. Determination of the node closest to the p ericenter of the initial orbit. Calculation of the time of passing this node. Op erations described by items 1.1 , 1.2 and 1.3 are performed fo r this purpose. 1.1. Calculation of the elements of the initial orbit: a, e, , i, , argumen t of latitude

u using state vector ( rE1 , v
1.2.
0

E1

)T

for the ti me t
0

E1

. : (68)

Calculation of tru e ano malies

and

for arguments of lati tude 0 and

=

,

=

27


Adjustment o f

0

and

to the interval:

,

. Calculation of the argument of latitude

for the node which is closer to the pericenter using the formula:

u1i =

0, ,

if

0

<

,

otherwise.

(69)

1.3. Calculation of the time for which the argu ment of l atitude is u1i . For this purpose we calculate t1i the time of the trans fer fro m the point with argument of latitude u1i to the point with argu ment of latitude u. Th e time fo r which t he argu ment of latitude is equal to u1i is calculated by the formula: t1i = t1 t1i . 2. Calculation of the period of the initial orbit by the fo rmula:

T1 = 2

a3 µ

(70)

3. Preliminary calculation of the time interval of the transfer. Operations 3.1 ­ 3 .4 are p erformed for this purpose. 3.1. Calculation of the state v ector o f the initial orbit fo r the time t1i : r t1i , v t1i 3.2. Calculation of the orbital elements: a2 , e2 ,
2

( ( ) ( ))

T

.

, i2 ,

2

using the state vector

(

rE 2 , v

E2

)T

for the time t

E2

.

3.3. Calculation of the semi-major axis of the transfer o rbit: a2 + r ( t1i ) atrn = 2 3.4. Calculation of the time of transfer:
3 atrn

(71)

ttrn =

µ

(72)

4. Evaluation of the times of the 1st and the 2nd thrusts. Operations described in 4.1 -4.5. 4.1. Enu meration of the possible times of the 1st thrust: t1k = t1i + kT1 , k = 2, 1, 0, 1, 2, 3 4.2. det ermination of the time of the 2nd thrust: t2 k = t1k + Ttrn 4.3. Calculation of the state v ector r1k , v of the state vector r2 k , v

(73) (74)

(

2k

)

(

1k

)
2k

for the time t1k using initial orbit. Calculation using the target orbit.

for the time t

Det ermination of the elements of the orbit for the transfer from position r1k for the time

t1k to the position r2 k for the ti me t

Calculation of the state vector r1tk , v

(

r2tk , v

2tk

)

(

2k

using Lambert's p roblem procedure.

1tk

)

for the time t1k and the stat e vector

for the time t

2k

using the elements of the transfer orbit .

Calculation of the module of the resulting impulse:

28


mk = v

1k

v

1tk

+v

2k

v

2 tk

(75)

4.4. Selection of the k , fo r which we h ave the mini mum mk . Th e times corresponding to

mk will be denoted as t1m and t2 m . 5. Updating of the times of the 1st and the 2nd thrusts
5.1. Enu merative search for t1a within t1m 5.2. Enu merative search for t within t
2m
2a

ttrn 4 ttrn 4

, t1m + ,t
2m

ttrn 4 ttrn 4

with the step 5 min utes. with the step 5 min utes.

2a

+

5.3. For each pair o f the v alues of t1a and t

we perform the operations similar to those

of item 4.3. We cal culate the state vector for the initial orbit for the time t1a and the target one ­ for the time t2 a . Then solve the Lambert's problem and determine the transfer impulses. Then calculate the module of resulting impul se. 5.4. Select the pair of ti mes t1a and t2 a , fo r which the minimum of resulting impulse is attained and gen erate the result. ALGORITHM FOR EVALUATI ON OF THE ACC ELERATI ON AND THE TI MES OF THE SWITCH-ON AND CUTOFF OF THE LOW THRUST ENGI NE BY TWO I NITI AL CONDITI ONS This section describes the algorithm for d etermining the time interval of the operation of the low thrust engine and the generated acceleration provid ing the transfer bet ween the given orbits. Input data: , ( rE1, v
2L

(
t

t

E1

t1L , t
E2

)

E1

)T

­ ­

time and the state v ector o f the initial orbit for this time; interval of possible engine s witch-on; time and the state v ector o f the t arget orbit for this time; interval of possible engine cutoff; we assume that the end of the interv al of possible switch-on is does not exceed the b eginning of the cutoff interv al; maximu m possible value of the step.

(

, rE 2 , v
2R

(

t1R , t

)

E2

)T

­ ­ ­

hmax

Output information:

t t

I1 I2

­ ­ ­ ­

estimate of the s witch-on ti me; estimate of the cutoff time; sign of coordinate frame: 0 ­ inertial, 1 ­ orbital;. acceleration vector.

Frame

a = ax , a y , a
Designation

(

z

)

T

29


y S t = rS t , v S t x
T T T

() ( () ()) t t t () = (r () , v ())

T

­ ­

state vector o f initial orbit for the time t ; state vector o f targ et orbit for the time t ;

T

Basic relationships The time within the interval of possible switch-on t1L , t time of cutoff fro m the interv al t1R , t

(

2R

)

(

2L

)

, will be denoted as t
R,i

L, j

, and the

-- interval of possible cutoff - as t
L, j

.

For the set times of s witch-on and cutoff of lo w thrust engines: t

,t

R,i

we d etermine the ac-

celeration vector under the condition of minimum weighted residual between the state vectors:

xi = x

T

( t R,i )

and
t
R ,i

yi +
t
L, j

(

t

R,i

,

) B( )

d

a

(76)

where

x y

i i

­ ­

state vector o f the sp acecraft in the target orbit for the ti me t

R,i

;
T

state vector o f the sp acecraft in the initial orbit for the t ime t Ri , y i = y

( t R, i )
0 E3

;

matrix 6x3, describing the direction of the thrust; if the direction of the thrust vector is permanent in the inertial coordinate frame, the matrix B

B

()

()

=

;

­

if the direction of the thrust is permanent in the orbital coordinate frame,

B E
3

()

=

0 er , e n , e

;
b

­
b

unit matrix of the 3rd order 3; unit vectors of the directions of the axes of RNB coo rdinate frame; the sought accel eration vector; matrix of partial derivatives of the state vector for the time t with respect to the state vector for the time of initial conditions along the i nitial orbit.

er , e n , e a

yt yt

( E1 )

()

­ ­ ­

(

t

R,i

,

)

=

yt

()
( E1 )
R,i

y

yt

yt

( E1 )

()

1

­ transition matrix

The six-dimensional diagonal matrix Wg is used as the weighting matrix in the search of the minimum. Diagonal elements of this matrix are defined by the values of a priori RMS errors of the position and velocity components: P and V :

30


1 Wg =
2 P

E3

0 1

3

(77)

0
Here 03 - zero matrix of the 3rd order.

3

2 V

E3

In case we use TLE the v alues of thes e RMS errors are: If we denote
t
R ,i

P

= 100 m/s,

V

= 1 m/s.

Q

j ,i

=
t
L, j

(

t

R,i

,B

)()

d

(78)

the considered weighted residual is repres ented in shap e of the follo wing quadratic form:

F= x

(

i

y

i

Q j ,ia

)

T

Wg x

(

i

y

i

Q j ,ia

)

(79)

For the fixed indexes i, j the vector a , for which the quadratic form (79) reach es its minimum, can be found using the formula:

a

min

= QT,i Wg Q j

(

j ,i

)
)
T

1

QT,i Wg x j

(

i

y

i

)

(80)

The respective minimum value of the qu adratic form is equal to:

Fmin = x

(

i

y

i

Q j ,ia

min

Wg x

(

i

y

i

Q j ,ia

min

)

(81)

The global minimu m o f the qu adratic form (79) can be found by scanning in i and j . For reduction of computation for the enumerative search we need the recursive fo rmulas connecting matrices Q j ,i +1 and Q j 1,i with the matrix Q i, j . Let us derive these formulas. Actually,
t
R ,i+1

Q
t

j , i+1

=
t
L, j

(

t

R , i+1

,

) ()
B =

d=

(

t

R , i+1 R , i t

,t

)

t

R ,i

t

( (
t

t

R,i

,

) B ( )d ) B ( )d

+
(82)

L, j

R , i+1

t

(

t

i+1

,

) B ( )d

(

t

R , i+1 , t R , i

)

R ,i+1

Q

j ,i

+
t
R ,i

R , i+1

,

i

31


t

R ,i

Q

j 1, i

=
t
L, j 1

(
t
L, j

t

R,i

,B

)()

t

L, j

d=
t
L, j 1

(

t

R,i

,B

)()

t

R ,i

d+
t
L, j

(

t

R,i

,B

)()

d=
(83)

(
Using
t
R , i+1

t

R,i

,t

L, j

)

t

(

t

L, j

,B

)()

d +Q
for

j ,i

L, j 1

the

trapezoid

formulas
t
L, j

app roximate

calculation

of

the

integrals

t

(

t

R , i +1

,

) ()
B

d

and
t
L, j 1

(
+ t

t

L, j

,B
t

)()
R,i

d yields:

R, i

Q

j , i +1

(
1 2

t

R , i +1 R , i

,t

)

Q

R , i +1

j ,i

2
t
Lj

Bt

(

R , i +1

)+ (
Lj

t

R , i +1 , t R , i

) B (t )
R,i

(84) (85)

Q

j 1, i

(

t Ri , t

Lj

) ( B (t ) + (
Lj

,t

Lj 1

) B ( t )) (
(

t

Lj

t

Lj 1

)

+Q

j ,i

Sequence of the operations
1. Select the numb er o f fragmentations N of the interval t1R , t under the condition: ble engine cutoff:
R (86) N1 2. Within the interval t1R , t2 R using parameters of the target orbit with perman ent step hR calculate the table of state v ectors:
2R

)

of possible engine cutoff

t

2R

t1

R

N1

< hmax and calculate the v alue of the step for the interval of possit
2R

hR =

t1

(

)

{

xi = x

T

( t R,i )

, t R, i = t1R + ( i 1) hR , i = 1, ...N

}

(87)

3. For each time t R ,i , i = 1, ..., N using paramet ers of the initial orbit we calculate the matrix o f the derivatives of the co mponents of the state v ector for the time t R ,i with respect to the co mponents of the state vector for the initial time t E1 :

Zt

()
R, i

=

yt

()
( E)
R, i

y t1

(88)

4. Then select the numb er of fragmentations M of the interv al of possible s witch-on of the en t t gine t1L , t2 L under the condition that the step o f the frag mentat ion hL = 2 L 1L does not M1 exceed hmax .

(

)

32


5. Within the interval of possible engine switch-on (t1L , t2 L ) using the parameters of the initial orbit with perman ent step hL we calculate the table co mprising the time

t

L, j

= t1L + ( j 1) hL , j = 1, ..., M , state v ector of the initial orbit for this time and the matrix of
E1

partial derivatives of the co mponents of the state vector for the initial time t components of the state v ector for the time t
1
L, j

with respect to the

t

L, j

, y j, Z

j

=

yt

()
( 1) E
L, j

1

yt

=

yt yt

()
L, j

( E1 )

Further operation described b y items 6-12 are perfor med for two variants of the direction of the thrust vector: perman ent in the inertial coordinate frame and per man ent in the orbital coordinate frame. 6. Calculate the matrix Q
M ,1

by performing the following operations. If t
M ,L

M ,L

= t1R , then matrix

Q

M ,1

= 0 . If t

M ,L

< t1R , the interval t

, t1

R

is frag mented into several p arts in a way that

the length of no p art ex ceeds hmax (see the pro cedure for selecting the step in items 1 or 3). The
t1
R

matrix Q

M ,1

=
t
ML

(

t1 R ,

) B ( )d

is calculated using trapeziu m fo rmula for each part.

7. Det ermine the v alues of the follo wing variables: ­ current minimum value of the quadratic form, Fc ­ acceleration vector corresponding to the current minimum value of the quadratic form, a
c

ic , jc

­

indices, corresponding to the current minimum value o f the quadratic form.

Initial value Fc =+ 8. We organize t wo cycles. The external cycle for j from M until 1 . Th e internal cycle for i from 1 until N . 9. Prior to begin o f the internal cy cle we set the variable matrix Q c = Q j 1 . 10. Internal cy cle. Using the matrix Q c , by Eqs (80) and (81) we find the acceleration vector a and the value of functional F . If F < Fc , new v alues are assigned to the variables Fc , a c , ic , jc . Using Eq. (84) and the current matrix Q c = Q
j ,i

we cal culate the new matrix Q c :

Qc = Q

(89) 1 B t R , i +1 + t R , i +1 , t R , i B t R , i t R , i +1 t R , i R , i +1 Ri j ,i 2 11. After the co mpletion of the intern al cycle, using the recurrent formula (4.10) we calculate the new matrix Q c :

(

t

,t

)

Q

+

((

j , i +1

=

)(

) ( )) (

)

1 2

(

t R1 , t

L, j

)( ( ) (
Bt
L, j

Qc = Q t
L, j

j 1,1

=
L, j

+

,t

L, j 1

) ( )) (
Bt

t

L, j

t

L, j 1

)

+Q

(90)
j ,1

33


12. We co mpare the values of Fc , cal culated for the t wo variants of the direction of the thrust vector. Th en select the minimum value and gen erate th e result using a c , ic , jc . REFERENCES

1 2 3

M.F. Subbotin, A Course of Celestial Mechanics, (in Ru ssian), Moscow, Len ingrad, vol. 1, 1933 P.E. E liasb erg, Introduction to th e Theo ry of Motion of Earth S atellites, "N auka", Mo scow , 1965 (in Ru ssian ).

D. E. Okhotsimsky and Yu. G. Sikharulidze. Foundations of th e Space Flight mech anics. (in Russian), "N auka", Moscow, 1990.
4 5

D.A. Vallado. Fundam entals of Astrodynamics and Applications. T hird Edition. Microco sm Press and Springer, 2007.

L. Euler. Determnatio orbitae com etae qui m ense martio huius anni 1742 potissimum fuit observatu s, Miscellanea Berolinensia 7, 1743.
6 7

J.H . L ambert. In signiores O rbita Com etarum Proprietates. Augusta Vindelico rum, Aug sburg, 1761.

A. Albouy. Lecons dur le p robleme d es deux corps. Collected p apers: Kepler's problem. Collisions. Regularization. Izh evsk, Institu te of computer research, 2006.
P.S. Laplace. Sur le prin cipe de la gravitation inierselle et sur les inegalites secularies des planetes qui en dependet, 1773j, oeuvres completes, G authier-Villars, Paris, tome 8, p 214.
9 8

J.C . Ad ams. On a S imple Proof of L ambert's Theorem, Messenger of Math . 7(1877) pp. 97-100. O. D ziob ek. Math em atical Theories of Planetary Motions. 1888. E.J. Routh. A Treatise on Dynam ics o f a Particle, C ambridge U , Press 1898. H.C. Plumm er. An Introductory Treatise on Dynam ical Astronomy, London, 1918. R.H. Battin. A stronomical Guidan ce, McGraw-Hill, New York, 1964.

10 11 12 13 14

V. Agapov, I. G alustov, S. K amen sky, A. Tuchin. Analysis o f Lau nch Profiles with Change o f In clin ation Based On Data On P rim ary Determin ed Orb its. P rep rint Inst. Appl. Mathem . Russia A cademy o f S ciences, 2007, 40.

34