Forming and keeping fast fly-around under constant thrust

Available online at www.sciencedirect.com

Advances in Space Research 48 (2011) 1421–1431 www.elsevier.com/locate/asr

Forming and keeping fast fly-around under constant thrust Yongqiang Qi ⇑, Yingmin Jia The Seventh Research Division, Beihang University, 37 Xueyuan Rd, Haidian District, Beijing 100191, China Received 9 March 2011; received in revised form 21 May 2011; accepted 15 June 2011 Available online 22 June 2011

Abstract A new switching control algorithm under constant thrust is designed for the chaser fast flying around the target spacecraft along a specified fly-around trajectory. The switching control laws are obtained based on the acceleration sequences and the on time of thrusters which can be computed by the time series analysis method. The perturbations and fuel consumptions are addressed during the computation of the on time of thrusters. Furthermore, the relative position parameters of the target spacecraft are obtained by using the vision measurement and the target fly-around positions are calculated through the isochronous interpolation method. The change of the relative position and the relative velocity of the chaser during the constant thrust fast fly-around are presented through simulation example. It is proved that, with the switching control laws, the chaser will fast fly around the target spacecraft along the specified fly-around trajectory. Crown copyright Ó 2011 Published by Elsevier Ltd. on behalf of COSPAR. All rights reserved. Keywords: Fast fly-around; Constant thrust; Switching control laws

1. Introduction Fast fly around formation (FFA) means the chaser can fly around the target spacecraft more than once in an orbital period under control. FFA can be widely used in space missions such as space rendezvous and docking (RVD) (Wang et al., 2009). In natural situations, the formation period is equal with the target spacecraft’s orbital period which is usually more than 1.5 h. Actually, some special missions such as space station maintenance need the chaser fly around the target spacecraft in a very short time (Enrico, 2008). The issues of controlling the FFA formation were discussed in previous research. Masutani et al. (2001) presented the bi-elliptic fly around maneuver by impulsive thrust control, and proposed an optimal feedback control scheme for the thrust to maintain this trajectory in the presence of disturbances. Straight et al. (2004) investigated a guidance scheme for FFA maneuvers involving impulsive burns at specific points within a specified toroidal region ⇑ Corresponding author.

E-mail addresses: [email protected] (Y. Qi), [email protected] (Y. Jia).

centered on the circular-orbiting reference satellite, and presents two analytical methods for determining the magnitude and direction of the impulses. Luo et al. (2006, 2007) developed the trajectory design and guidance algorithms for entry-into circumnavigation, fast circumnavigation and exit circumnavigation based on multiple impulse maneuver method. Wang et al. (2007) studied the initialization condition for chaser flying around target in eccentric orbits, and designed LQR control laws to keep chaser in scheduled cycle. Zhu (2009) described five fly around patterns named natural elliptical, natural spiral, single-impulse tear-drop, multi-circular, and multi-impulse gymkhana. There are two ways to realize fast fly around formation, one is impulsive control and the other is continuous variable thrust control (Kojima, 2003; Zhou et al., 2006; Zachary, 2008). Most of the researches about FFA mentioned above are based on impulse thrusts (Wu et al., 2006; Yasuhiro et al., 2001; Zhou et al., 2006). In actual practice, however, the thrusters of the spacecraft are of constant thrust, therefore, fast fly around cannot normally be considered as continuous thrust control or impulsive control. In this paper, we build the dynamic model of fast space circle formation using constant thrust control, and derive analytical expressions to compute minimum fuel expenditure during a

0273-1177/$36.00 Crown copyright Ó 2011 Published by Elsevier Ltd. on behalf of COSPAR. All rights reserved. doi:10.1016/j.asr.2011.06.013

1422

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431

formation period. The control strategy presented in this paper can not only be used in fast space circle formation on circular orbits, but also be appropriate to any kind of fast fly around formations on eccentric orbits. The purpose of this paper is to realize fast fly around formation along the specified fly-around trajectory based on constant thrust control. The specified fly-around trajectory is given according to the initial relative position and the terminal target position of the chaser, at the same time the size and the direction of the initial relative velocity and the acceleration are also taken into account (Richards et al., 2002; Schouwenaars et al., 2004). First of all, the relative position parameters of the spacecraft are obtained by using vision measurement. Next, the definition of the specified fly-around trajectory is presented. Then, to ensure the chaser flying along the specified trajectory, the isochronous interpolation method is used by dividing the entire process of fly-around into equal time intervals. At last, the on time of thrusters can be respectively computed by using the time series analysis method. Moveover, the switch control laws for forming and keeping fast fly-around are designed based on the on time series of thrusters corresponding accelerations. The simulation results show that the chaser flying along the specified fly-around trajectory under the constant thrust. 2. Relative position determination algorithm based on vision measurement The purpose of vision measurement is to get the relative position parameters of the target spacecraft. Two cameras C1 and C2 which have the same focal lengths and the internal parameters are installed in the chaser. Camera coordinate system is established as shown in Fig. 1. Oc1 and Oc2 are the optical centers of the cameras C1 and C2, respectively. X c1 and X c2 are the optical axes of the cameras C1 and C2, respectively. X c1 is parallel to X c2 , and Y c1 coincides with Y c2 , where the optical parallax between Oc1 and Oc2 is L. P is any point on the target spacecraft, P1 and P2 are the projection points of the P on the image planes of the cameras C1 and C2, respectively. Therefore, P is the intersection of Oc1 P and Oc2 P .

Further assume that the cameras have been calibrated and their projection matrices are M1 and M2 (Yamanaka and Ankersen, 2002). The following results can be obtained by coordinate transformation 0 1 0 1 0 1 0 1 X X x1 x 2 BY C BY C B C B C B C B C X c1 @ y 1 A ¼ M 1 B C; X c2 @ y 2 A ¼ M 2 B C; ð1Þ @ZA @ZA 1 1 1 1 0 i 1 m11 mi12 mi13 mi14 B i C M i ¼ @ m21 mi22 mi23 mi24 A; i ¼ 1; 2; ð2Þ mi31

Oc1

OC

P1 Z c1

Yc2

Oc2 Zc2

mi34

ðx1 m131  m111 ÞX þ ðx1 m132  m112 ÞY þ ðx1 m133  m113 ÞZ ¼ m114  x1 m134 ;

ð3Þ

ðy 1 m131  m121 ÞX þ ðy 1 m132  m122 ÞY þ ðy 1 m133  m123 ÞZ ¼ m124  y 1 m134 ;

ð4Þ

ðx2 m231  m211 ÞX þ ðx2 m232  m212 ÞY þ ðx2 m233  m213 ÞZ ¼ m214  x2 m234 ;

ð5Þ

ðy 2 m231  m221 ÞX þ ðy 2 m232  m222 ÞY þ ðy 2 m233  m223 ÞZ ¼ m224  y 2 m234 :

ð6Þ

The linear equations can be rewritten in the form of matrix equation (Philip and Ananthasayanam, 2003) T

AðX ; Y ; ZÞ ¼ B; 0 x1 m131  m111 B y m1  m1 B 21 A ¼ B 1 231 @ x2 m31  m211 y 2 m231  m221 0 1 m14  x1 m134 B m1  y m1 B 24 1 34 B 2 @ m14  x2 m234 m224  y 2 m234

ð7Þ x1 m132  m112 y 1 m132  m122 x2 m232  m212 y 2 m232  m222 1 C C C: A

1 x1 m133  m113 y 1 m133  m123 C C C; 2 2 A x2 m33  m13 y 2 m233  m223

ð8Þ

The least-square solution of P in the reference coordinate system as follow (Kukush and Huffel, 2004),

Xc2 Yc1

mi33

where (x1, y1, 1)T and (x2, y2, 1)T are homogeneous coordinates of P1 and P2 in the cameras C1 and C2, respectively. (X, Y, Z, 1)T is the homogeneous coordinate of P in the reference coordinate system, that is to say, X, Y and Z are the 3D position coordinates of the real point P. Therefore, X c1 and X c2 can be removed by combining (1) and (2),

B¼

Xc1

mi32

P2

L

Fig. 1. Model of binocular vision.

T

1

ðX ; Y ; ZÞ ¼ ðAT AÞ AT B:

ð9Þ

Assume that the virtual optical center of the cameras C1 and C2 is the midpoint of the line segment OC1 OC2 and defined as OC, and the virtual optical axis is parallel to Oci X ci ði ¼ 1; 2Þ, the other two axes are coincide with

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431

Oci Y ci ði ¼ 1; 2Þ and parallel to Oci Z ci ði ¼ 1; 2Þ, respectively. According to the definition of the coordinate system (Geller, 2006), fL x1  x0 ; Y ¼ L; ðx2  x1 Þdy ðx2  x1 Þ ðy  y 0 Þdz Z¼ 1 L; ðx2  x1 Þdy X ¼

ð10Þ

where f is the focal length of the camera, (x0, y0) is the projection point of the P on the image plane from the virtual optical center OC. dy and dz are the physical sizes in the y-axis and the z-axis represented by one pixel on the image plane. Therefore, the following results can be obtained according to vision measurement errors @X @X dX ¼ dx1 þ dx2 @x1 @x2 fL fL dx1  dx2 ; ¼ 2 ðx2  x1 Þ dy ðx2  x1 Þ2 dy @Y @Y dx1 þ dx2 @x1 @x2 ðx2  x0 ÞL ðx1  x0 ÞL ¼ dx1  dx2 ; 2 2 ðx2  x1 Þ ðx2  x1 Þ

ð11Þ

dY ¼

@Z @Z @Z dx1 þ dx2 þ dy @x1 @x2 @x2 1 ðy  y 0 ÞdzL ¼ 1 ðdx1  dx2  dy 1 Þ; 2 ðx2  x1 Þ dy

ð12Þ

dZ ¼

ð13Þ

where the independent parameters d x1, d x2, d y1, d y2 are the pixel position measurement errors of the cameras C1 and C2, so it can be assumed that d x1 = d x2 = d y1 = d y2, therefore, they have the same variance D(x1) = D(x2) = D(y1) = D(y2) = D(x). Next, the position errors can be obtained by calculating the variances on both sides of the Eqs. (11)–(13). 

@X @x1 "

DX ¼ ¼2

@Y @x1

2



ðx2  x1 Þ  "

@Z @x1

4

@Y @x2

 Dðx1 Þ þ

ðx2  x1 Þ

ð14Þ

2 Dðx2 Þ ð15Þ

2  2 @Z @Z Dðx2 Þ þ Dðy 1 Þ @x2 @x2 # 2

4

Dðx2 Þ

i 2 2 ðx2  x0 Þ þ ðx1  x0 Þ ;

2

3L2

2

DðxÞ;

Dðx1 Þ þ

L2 DðxÞ h

DZ ¼

@X @x2

#2

fL 2



¼

 Dðx1 Þ þ

ðx2  x1 Þ dy

DY ¼ ¼

2

ðy 1  y 0 Þ DðxÞ:

ð16Þ

1423

It is clear that position errors DX, DY and DZ are inversely proportional to optical parallax. In order to ensure the target spacecraft always within the range of vision, the focal lengths of the cameras should be adjusted dynamically to change the size of the range of vision. 3. Fast fly around formation design based on constant thrust The process of forming and keeping fast fly-around is divided into two stages: the coplanar adjustment stage and the fast fly-around stage. 3.1. The coplanar adjustment stage After the relative position parameters of the target spacecraft are obtained by using vision measurement, then the position parameters of the chaser relative to the target spacecraft can be calculated through coordinate transformation. Therefore, the relative motion coordinate system can be established as follows: the target spacecraft is assumed as a rigid body on a circular orbit; OT is selected as the origin of coordinate; x-axis is along the direction opposite to the target’s velocity vector, y-axis is along the radius direction from the earth’s central mass to the target’s central mass and z-axis completes the right-handed frame. Therefore, the relative motion can be described by Clohessy–Wiltshire equations. €x  2xy_ ¼ ax þ apx ;

ð17Þ

€y þ 2x_x  3x2 y ¼ ay þ apy ;

ð18Þ

€z þ x2 z ¼ az þ apz ;

ð19Þ

where x represents the mean angular velocity of the target spacecraft. ax, ay and az represent the thrust accelerations of the chaser in three axes, respectively. It is clear that the relative motion can be decomposed into the motions in the xoy-plane and the z-axis. Due to the deviations of the earth’s gravitational potential from a sphere and the atmospheric drag and solar pressure as well as thruster plume pressure etc, the nonlinear perturbations in the Clohessy–Wiltshire equations should be taken into consideration. The most significant effect on near Earth orbits is caused by the second harmonic of the Earth potential, which represents the Earth’s oblateness. The coefficient J2 is more than two orders of magnitude larger than all other ones, hence the name J2-effect for the most significant gravitational disturbance (Fehse, 2003). The effect of the oblateness of the Earth on the orbits results in a motion of the line of nodes called the “regression of nodes”. Wertz and Larson (1991) given the following formulas for the rate of change of the line of nodes  2 _XJ ¼ 1:5xJ 2 RE ðcos iÞð1  e2 Þ2 ; ð20Þ 2 a

1424

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431

where RE is the Earth radius at the equator, a is the semimajor axis in kilometres, e is the eccentricity and i is the inclination. Solar radiation produces a force on a spacecraft in the Sun-spacecraft direction. As the forces due to solar pressure have in- and out-of-plane components, depending on the Sun direction w.r.t the orbital plane, the solar pressure will have some effect on all orbital parameters, with the most important ones being on eccentricity and on inclination. Depending on orbital height and Sun direction, the force will be intermittent, i.e. the force will be zero when the satellite is in the shadow of the Earth. Solar pressure is the most prominent disturbance for rendezvous trajectories in geosynchronous orbits, where drag is practically zero but where solar pressure, in combination with a difference in the ballistic coefficients of chaser and target vehicles, can lead to different accelerations of the two. Because the Sun-spacecraft direction varies along the orbit and with the year, the actual effects of solar pressure need to be calculated separately for each case. In a more detailed analysis, and for verification purposes, the individual surfaces of the satellite body have to be taken into account concerning their direction w.r.t. the Sun-spacecraft vector and their reflection properties, i.e. whether absorptive, diffuse reflective or specular reflective (Fehse, 2003). In a 356 km altitude orbit, the acceleration due to solar pressure and the resulting effect on the rendezvous trajectories is about two to three orders of magnitude lower than that of atmospheric drag. Suppose that apx ; apy and apz represent the sum of the perturbation and nonlinear factors in three axes, respectively. And apx ; apy and apz satisfy the following conditions: apx  gx ;

apy  gy ;

apz  gz

ð21Þ

Although gx, gy, gz are variables, but they vary little during the fly-around manoeuvres and can be approximately seen as constants. Therefore, ax, ay, az are correlative only with the size of the thrusts and the mass of the chaser: ai ¼ 

Fi ; m0  ðm_ x tx þ m_ y ty þ m_ z tz Þ

i 2 fx; y; zg;

ð22Þ

where Fx, Fy, Fz represent the vacuum thrusts of the chaser in three axes, respectively. m0 represents the initial mass of the chaser at the beginning of the coplanar adjustment stage. tx, ty, tz represent the on time of thrusters and m_ x ; m_ y ; m_ z represent the mass-flow-rate of propellant of the chaser’s thrusters in three axes, respectively. Suppose that in the process of the coplanar adjustment the thrust accelerations in the x-axis and the y-axis are zero and there is no possibility of collision between the chaser and the target spacecraft. Therefore, only the relative motion in the z-axis need to be taken into consideration. The form of the polynomial approximate expression of Eq. (22) is: az ¼

n X Fz ¼ Fz kzj ðT c Þj ; m0  m_ z tz j¼0

where kzj, j 2 {1, 2, . . . , n} are the polynomial coefficients and determined by the number of thrusters fired and their working times. k0 = 1/m0 represents the inverse of initial mass of the chaser at the beginning of the coplanar adjustment stage. Tc is the time of the coplanar adjustment stage and tz 6 Tc. Then the original system can be transformed into the form: €z þ x2 z ¼ gz þ F z

n X

j

kzj ðT c Þ :

ð24Þ

j¼0

The coplanar adjustment is mainly determined by Eq. (24) in which z stands for the extent of non coplanarity. The task of the coplanar stage is to make the chaser coplanar with target spacecraft as shown in Fig. 2. The key of coplanar adjustment is to calculate the on time of thruster in the z-axis. The coplanar adjustment can also divided into three stages: the accelerating stage, the zero-thrust stage and the decelerating stage. Suppose that taz is the time of accelerating stage, tbz is the time of zero-thrust stage and tcz is the time of decelerating stage. tz ¼ taz þ tcz ;

ð25Þ

T c ¼ taz þ tbz þ tcz :

ð26Þ

The target relative position and the target relative velocity of the coplanar adjustment both are zero. Then the following two results can be obtained according to Eq. (24)   Fz 2F z 0 ¼ z0  2 ðkz0 þ gz Þ þ 4 kz2 cos xT c x x   vz0 F z Fz  kz1 sin xT c þ 2 ðkz0 þ gz Þ þ x x3 x 2F z Fz Fz 2 ð27Þ  4 kz2 þ 2 kz1 T c þ 2 kz2 ðT c Þ ; x x x   Fz 0 ¼ vz0  2 kz1 cos xT c x   Fz 2F z Fz  xz0  ðkz0 þ gz Þ þ 3 kz2 sin xT c þ 2 kz1 x x x Fz þ 2 2 kz2 T c ; ð28Þ x where z0 and vz0 are the initial relative position and the initial relative velocity of the chaser in the z-axis at the begin-

ð23Þ Fig. 2. The chaser’s trajectory in the coplanar stage.

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431

ning of the coplanar adjustment. Then the on time of thruster in the z-axis can be calculated as follows: (1) When z0vz0 > 0, the time of three stages in the z-axis can be calculated by the following equations F z tcz1 ; jvz0 j ¼ m0  m_ z tcz1 taz tcz2 ¼ ; m0  m_ z ðtcz1 þ taz Þ m0  m_ z ðtcz1 þ tcz2 þ taz Þ jz0 j þ

ð29Þ ð30Þ

F z t2cz1 F z t2cz2 ¼ 2ðm0  m_ z tcz1 Þ 2½m0  m_ z ðtcz1 þ tcz2 þ taz Þ þ

F z t2az ðtaz þ 2tbz Þ ; 2½m0  m_ z ðtcz1 þ taz Þ

ð31Þ

where tcz1 and tcz2 are the time of two different decelerating stages in the z-axis. tcz1 is the initial stage of the coplanar adjustment and tcz2 is the final stage of the coplanar adjustment. Then the time of coplanar adjustment t can be obtained, T c ¼ tcz1 þ taz þ tbz þ tcz2 :

ð32Þ

(2) When z0vz0 6 0, the time of three stages in the z-axis can be calculated under two different conditions. If z0 and vz0 meet the following conditions jvz0 j ¼

F z tcz1 ; m0  m_ z tcz

jz0 j <

ðm0  m_ z tcz1 Þv2z0 ; 2F z

ð33Þ

then the time of three stages in the z-axis and the time of coplanar adjustment t can be calculated as follows taz tcz2 ¼ ; m0  m_ z ðtcz1 þ taz Þ m0  m_ z ðtcz1 þ tcz2 þ taz Þ

F z t2az ðtaz þ 2tbz Þ ; 2½m0  m_ z ðtcz1 þ taz Þ

T c ¼ tcz1 þ taz þ tbz þ tcz2 :

F z tcz1 ; m0  m_ z tcz

jz0 j P

ð35Þ ð36Þ

ðm0  m_ z tcz1 Þv2z0 2F z

ð37Þ

in this case, there are only the zero-thrust stage and the decelerating stage during coplanar adjustment. The time of the zero-thrust stage and the decelerating stage and the time of coplanar adjustment t can be calculated as follows jz0 j ¼

F z t2cz þ jvz0 jtbz ; 2ðm0  m_ z tcz Þ

T c ¼ tbz þ tcz :

The fast fly-around stage consists of two steps. The first step is to make the chaser enter the relative fast fly-around orbit. The second step is to keep the fast fly-around. The manoeuvres of forming and keeping fast fly-around are carried out in the xoy-plane. Therefore, the form of the polynomial approximate expression of Eq. (22) is: ax ¼

N X Fx j ¼ Fx kj ðT e Þ ; ðme  m_ x tx  m_ y ty Þ j¼0

ð40Þ

ay ¼

N X Fy ¼ Fy kj ðT e Þj ; ðme  m_ x tx  m_ y ty Þ j¼0

ð41Þ

where me is the initial mass of the chaser at the beginning of the first step. kj, j 2 {1, 2, . . . , N} are the polynomial coefficients and determined by the number of thrusters fired and their on time of thrusters during the first step. k0 = 1/me represents the inverse of initial mass of the chaser at the beginning of the first step. Te is the time the chaser takes to enter the relative fast fly-around orbit and tx, ty are the on times of thrusters in the x-axis and y-axis, respectively. Therefore, ti 6 Te, i 2 {x, y} and the Eqs. (17) and (18) can be transformed into the form: N X

j

kj ðT e Þ ;

ð42Þ

j¼0 N X

j

kj ðT e Þ :

ð43Þ

j¼0

If z0 and vz0 meet the following conditions jvz0 j ¼

3.2. The fast fly-around stage

€y þ 2x_x  3x2 y ¼ gy þ F y

F z t2cz1 F z t2cz2  jz0 j ¼ 2ðm0  m_ z tcz1 Þ 2½m0  m_ z ðtcz1 þ tcz2 þ taz Þ þ

different in different stages (the accelerating stage, the zero-thrust stage and the decelerating stage).

€x  2xy_ ¼ gx þ F x ð34Þ

1425

ð38Þ ð39Þ

After the time of coplanar adjustment Tc and the time of three stages taz, tbz, tcz are obtained, the polynomial coefficients kz0, kz1, kz2 can be calculated according to Eqs. (27) and (28). And kz0, kz1, kz2 are piecewise functions and are

The process of the chaser entering the relative fast flyaround orbit can also divided into three stages: the accelerating stage, the zero-thrust stage and the decelerating stage. Suppose that tax and tay are the time of accelerating stage, tbx and tby are the time of zero-thrust stage and tcx and tcy are the time of decelerating stage in the x-axis and in the yaxis, respectively. tx ¼ tax þ tcx ; ty ¼ tay þ tcy ; T e ¼ tx þ tbx ¼ ty þ tby :

ð44Þ ð45Þ

Suppose that the initial relative position of the chaser is (x0, y0) and the initial relative velocity of the chaser is (vx0, vy0) at the beginning of the chaser to enter the fast fly-around orbit. The target relative position in the fast fly-around orbit is (xT, yT) and the target relative velocity is (vxT, vyT). It is show that kj, j 2 {1, 2, . . . , N} are closely related to the mean angular velocity of the target spacecraft x. Due to x is small, kj is very smaller than x when j P 3 and have little impact on the polynomial approximate of Eqs. (42) and (43). Therefore, in the derivation process of

1426

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431

the formulas (46)–(53), kj(j P 3) can be ignored. Therefore, the following results can be obtained according to Eqs. (42) and (43), xT ¼ 2y 0 sin xT e  2ðvy0  b1 Þ cos xT e þ

4 X j¼2

If (vx0, vy0) and (vxT, vyT) meet the following condition     vxT  vx0  vyT  vy0   < ;  F   F  x

4 X 1 1 j j kj2 ðT e Þ 2x bj1 ðT e Þ þ Fx j ðj  1Þj j¼3

then the time of three stages in the x-axis and in the y-axis can be calculated by Eqs. (58)–(61).

þ x0 þ 2vy0  2b1 þ ðvx0  2y 0 xÞT e 1 þ ðk0 þ gx ÞF x T 2e : 2

vxT ¼ vx0 þ

F x tax ; m0  ðm_ x þ m_ y Þtax

ð58Þ

vyT ¼ vy0 þ

F y tax F y ðtay  tax Þ þ ; m0  ðm_ x þ m_ y Þtax m0  m_ x tax  m_ y tay

ð59Þ

ð46Þ

vxT ¼ 2xðvy0  b1 Þ sin xT e þ 2xy o cos xT e þ

3 X j¼1

F y t2ax ; jxT  x0 j ¼ vxo tax þ vxT tbx þ  2 m0  ðm_ x þ m_ y Þtax

3 X 1 j F x kj1 ðT e Þ 2xbj ðT e Þ þ j j¼2 j

þ vx0  2xy 0 þ ðk0 þ gx ÞF x T e : y T ¼ y 0 cosxT e þ

ð47Þ

4 X ðvy0  b1 Þ j sin xT e þ bj ðT e Þ ; x j¼1

vyT ¼ ðvy0  b1 Þ cos xT e  xy 0 sin xT e þ

4 X

jbj ðT e Þ

j1

; ð49Þ

b1 ¼ 

ð60Þ

F y tax ð2tay  tax Þ  jy T  y 0 j ¼ vyo tay þ  2 m0  ðm_ x þ m_ y Þtax 2

F y ðtay  tax Þ : þ  2 m0  m_ x tax  m_ y tay

ð48Þ

j¼1

2F x Fy 4F x ðk0 þ gx Þ þ 2 k1 þ 3 k2 ; x x x Fx Fy 2F x k2 ; b2 ¼  k1 þ 2 k2 ; b3 ¼  x x 3x

ð57Þ

y

ð61Þ

(2) When jvxTj P jvx0j and jvyTj < jvy0j, the time of three stages in the x-axis and in the y-axis can be calculated according to two different conditions.

ð50Þ F x k2 ; ð51Þ b4 ¼ 2x

where k0, k1, k2 are the above-mentioned polynomial coefficients and b1, . . . , b4 are expressed by k0, k1, k2. Then the on time of thrusters in the x-axis and the y-axis can be calculated as follows: (1) When jvxTj P jvx0j and jvyTj P jvy0j, the time of three stages in the x-axis and in the y-axis can be calculated according to two different conditions. If (vx0, vy0) and (vxT, vyT) meet the following condition     vxT  vx0      P vyT  vy0 ; ð52Þ  F   F  x y

If (vx0, vy0) and (vxT, vyT) meet the following condition     vxT  vx0      P vyT  vy0 ;  F   F  x

then the time of three stages in the x-axis and in the y-axis can be calculated by Eqs. (63)–(66). vxT ¼ vx0 þ

F x tcy F x ðtax  tcy Þ þ ; m0  ðm_ x þ m_ y Þtcy m0  m_ x tax  m_ y tcy

ð63Þ

vyT ¼ vy0 þ

F y tcy ; m0  ðm_ x þ m_ y Þtcy

ð64Þ

F x tay ð2tax  tcy Þ  jxT  x0 j ¼ vxo tax þ  2 m0  ðm_ x þ m_ y Þtcy

then the time of three stages in the x-axis and in the y-axis can be calculated by Eqs. (53)–(56). F x tay F x ðtax  tay Þ þ ; m0  ðm_ x þ m_ y Þtay m0  m_ x tax  m_ y tay F y tay vyT ¼ vy0 þ ; m0  ðm_ x þ m_ y Þtay F x tay ð2tax  tay Þ  jxT  x0 j ¼ vxo tax þ  2 m0  ðm_ x þ m_ y Þtay

vxT ¼ vx0 þ

ð53Þ

ð62Þ

y

2

F x ðtax  tcy Þ ; þ  2 m0  m_ x tax  m_ y tcy

ð65Þ

F y t2cy : jy T  y 0 j ¼ vyo tcy þ vyT tby þ  2 m0  ðm_ x þ m_ y Þtcy

ð66Þ

ð54Þ If (vx0, vy0) and (vxT, vyT) meet the following condition

F x ðtax  tay Þ2 ; þ  2 m0  m_ x tax  m_ y tay

ð55Þ

F y t2ay : jy T  y 0 j ¼ vyo tay þ vyT tby þ  2 m0  ðm_ x þ m_ y Þtay

ð56Þ

    vxT  vx0  vyT  vy0   < ;  F   F  x

ð67Þ

y

then the time of three stages in the x-axis and in the y-axis can be calculated by Eqs. (68)–(71).

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431

F x tax ; m0  ðm_ x þ m_ y Þtax F y tax F y ðtcy  tax Þ ¼ vy0 þ þ ; m0  ðm_ x þ m_ y Þtax m0  m_ x tax  m_ y tcy

vxT ¼ vx0 þ

ð68Þ

vyT

ð69Þ

F y t2ax ; jxT  x0 j ¼ vxo tax þ vxT tbx þ  2 m0  ðm_ x þ m_ y Þtax F y tax ð2tcy  tax Þ  jy T  y 0 j ¼ vyo tcy þ  2 m0  ðm_ x þ m_ y Þtax

ð70Þ

2

F y ðtcy  tax Þ : þ  2 m0  m_ x tax  m_ y tcy

ð71Þ

N axi X

tðaxiÞj þ

N cxi X

j¼1

tyi ¼

N ayi X

N cyi X

j¼1

T fi ¼

tðcxiÞj ;

ð72Þ

tðcyiÞj ;

ð73Þ

j¼1

tðayiÞj þ

N axi X j¼1

j¼1

tðaxiÞj þ

N cxi X j¼1

tðcxiÞj þ

N bxi X

tðbxiÞj ;

where mi is the initial mass of the chaser at the beginning of the ith thrust arc. where ki0, . . . , ki2 are polynomial coefficients in the ith thrust arc and ki0 ¼ m1i represents the inverse of initial mass the chaser at the beginning of the ith thrust arc. Suppose that the initial relative position of the chaser is (xi, yi) and the initial relative velocity of the chaser is (vxi, vyi) at the beginning of the ith thrust arc. The target position is (xi+1, yi+1) and the target velocity is (vx(i+1), vy(i+1)). Therefore, the following results can be obtained according to Eqs. (17), (18), (75), (76): xiþ1 ¼ 2y i sin xT fi  2ðvyi  bi1 Þ cos xT fi

Using the same method, the time of three stages in the xaxis and in the y-axis can be calculated for the case 3: jvxTj < jvx0j and jvyTj P jvy0j, and for the case 4: jvxTj < jvx0j and jvyTj < jvy0j. After the time the chaser takes to enter the relative fast fly-around orbit Te and the time of three stages (the time of accelerating stage tax and tay, the time of zero-thrust stage tbx and tby, and the time of decelerating stage tcx and tcy) are obtained, the polynomial coefficients k0, k1, k2 can be calculated according to Eqs. (45)–(51). k0, k1, k2 are piecewise functions and are different in different stages. The chaser begin to fly around the target spacecraft when it enters the fast fly-around orbit. The fast fly-around orbit is divided into N thrust arcs through the isochronous interpolation method and each thrust arc takes the same amount of time. Using the same method shown as above, the time of each thrust arc is divided into three types: the time of accelerating stages, the time of zero-thrust stages and the time of decelerating stages. Taking the working time of thruster in the x-axis as an illustration: suppose that Tfi is the fly time and txi is the on time of thruster in the ith thrust arc. Further assume that Naxi is the number of all the accelerating stages t(axi)j, j = 1, 2, . . . , Naxi. Ncxi is the number of all the decelerating stages t(cxi)j, j = 1, 2, . . . , Ncxi. Nbxi is the number of all the zero-thrust stages t(bxi)j, j = 1, 2, . . . , Nbxi in the ith thrust arc. Using the same method, the three types of stages can be defined in the y-axis in the ith thrust arc. txi ¼

1427

ð74Þ

þ

4 X j¼2

þ

4 X

1 2x biðj1Þ ðT fi Þj þ xi þ 2vyi  2b1 j Fx

j¼3

y iþ1

1 j kiðj2Þ ðT fi Þ þ ðvxi  2y i xÞT fi ðj  1Þj

1 þ ðki0 þ gx ÞF x T 2fi ; 2 4 X ðvyi  bi1 Þ j sin xT fi þ ¼ y i cosxT fi þ bij ðT fi Þ ; x j¼1

ð77Þ ð78Þ

vxðiþ1Þ ¼ 2xðvyi  bi1 Þ sin xT fi þ 2xy i cos xT fi þ

3 X

2xbij ðT fi Þj þ

j¼1

3 X 1 F x kiðj1Þ ðT fi Þj j j¼2

þ vxi  2xy i þ ðki0 þ gx ÞF x T fi ;

ð79Þ

vyðiþ1Þ ¼ ðvyi  bi1 Þ cos xT fi  xy i sin xT fi þ

4 X

jbij ðT fi Þ

j1

;

ð80Þ

j¼1

2F x Fy 4F x ðki0 þ gx Þ þ 2 ki1 þ 3 ki2 ; x x x Fx Fy bi2 ¼  ki1 þ 2 ki2 x x 2F x F x ki2 ki2 ; bi4 ¼ ; bi3 ¼  3x 2x

bi1 ¼ 

ð81Þ ð82Þ

where bi1, . . . , bi4 are coefficients in the ith thrust arc and can be expressed by ki0, . . . , ki2, ki0, . . . , ki2 can be obtained according to Eqs. (77)–(82). Then the thrust accelerations of the chaser ax and ay in the ith thrust arc in the x-axis and in the y-axis can be respectively calculated. In addition, in order to make the chaser fast flying along a circular flyaround orbit, the initial relative position (xi, yi) and the target position (xi+1, yi+1) of the chaser in the ith thrust arc should meet the following condition,  2  ðx þ y 2 Þ  ðx2 þ y 2 Þ 6 e; ð83Þ i i iþ1 iþ1

j¼1

ax ¼

2 X Fx j ¼ Fx kij ðT fi Þ ; ðmi  m_ x txi  m_ y tyi Þ j¼0

ð75Þ

ay ¼

2 X Fy ¼ Fy kij ðT fi Þj ; ðmi  m_ x txi  m_ y tyi Þ j¼0

ð76Þ

where e is the maximum permissible error of the chaser deflects from fly-around orbit. Further assume that Tfi is divided into Nxi and Nyi equal time intervals in the x-axis and in the y-axis, respectively. The equal time intervals are small enough so that the three types of stages (the accelerating stages, the zero-thrust

1428

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431

stages and the decelerating stages) can be equally divided. Then Tfi can be expressed as sequences of the equal time intervals in the x-axis and in the y-axis, respectively:

At last, the switch control laws in the ith thrust arc for fast fly-around can be given in the x-axis and in the y-axis, respectively:

fT xi1 ; T xi2 ; . . . ; T xiN xi g;

S x ¼ fðtxi1 ; axi1 Þ; . . . ðtxik ; axik Þ; . . . ðtxiN xi ; axiM xi Þg; n o S y ¼ ðtyi1 ; ayi1 Þ; . . . ðtyik ; ayik Þ; . . . ðtyiN yi ; ayiM yi Þ :

fT yi1 ; T yi2 ; . . . ; T yiN yi g:

ð84Þ

The on time of the thrusters can also be divided into Mxi, Myi(Mxi 6 Nxi, Myi 6 Nyi) equal time intervals in the ith arc. Therefore, the on time of the thrusters can also be expressed as sequences of the equal time intervals in the xaxis and the y-axis, respectively: ftxi1 ; txi2 ; . . . ; txiM xi g;

ftyi1 ; tyi2 ; . . . ; tyiM yi g

ð85Þ

because txi1 and tyi1 represent the on time of the thrusters in the first equal time interval and therefore they can be calculated as follows: txi1 2 f0; T xi1 g;

tyi1 2 f0; T yi1 g;

ð86Þ

suppose Txi1 = min{Txi1, Tyi1}, therefore, the appropriate value of txi1 and tyi1 can be selected by using cut and try method to make the following equation equal on both sides: 1 mi  m_ x txi1  m_ y tyi1

¼

2 X

kij ðT xi1 Þj :

ð87Þ

j¼0

In the same way, txi2 and tyi2 can be obtained 2 X 1 j

¼ kij ðT xi2 Þ : mi  m_ x ðtxi1 þ txi2 Þ  m_ y tyi1 þ tyi2 j¼0

ð88Þ

Therefore, ftxi1 ; txi2 ; . . . ; txiM xi g and ftyi1 ; tyi2 ; . . . ; tyiM yi g can be deduced by analogy. txik belongs to the zero-thrust stages if txik = 0, k 2 (1, 2, . . . , Mxi), but if txik – 0, it can be divided into the following three cases. ð1Þ jvxik j < jvxiðkþ1Þ j;

ð89Þ

vxik vxiðkþ1Þ P 0;

if vxik > 0, then txik belongs to the accelerating stages and Fx is in the positive direction of the x-axis. If vxik 6 0, then txik belongs to the accelerating stages and Fx is in the negative direction of the x-axis. ð2Þ jvxik j > jvxiðkþ1Þ j;

ð90Þ

vxik vxiðkþ1Þ P 0

if vxik > 0, then txik belongs to the decelerating stages and Fx is in the negative direction of the x-axis. If vxik 6 0, then txik belongs to the accelerating stages and Fx is in the positive direction of the x-axis. ð3Þ jvxik j > jvxiðkþ1Þ j;

ð91Þ

vxik vxiðkþ1Þ < 0;

if vxik > 0, then txik belongs to the accelerating stages and Fx is in the negative direction of the x-axis. If vxik 6 0, then txik belongs to the accelerating stages and Fx is in the positive direction of the x-axis. The thrust acceleration in the time interval txik in the x-axis is axik ¼ F x

2 X

j

kij ðtxik Þ ; ðtxik > 0Þ;

axik ¼ 0;

ðtxik ¼ 0Þ:

j¼0

ð92Þ

ð93Þ ð94Þ

When the magnitudes and the directions of thrust accelerations in the x-axis and the y-axis and the on time of the thrusters in each thrust arc are calculated, the switch control laws for the whole process of the fast fly-around can be obtained, then the goal of forming and keeping the fast fly-around motion can be achieved. 4. Numerical example Because the low earth orbits (LEO) carry orbit levels from 200–1200 km over earth’s surface. Therefore, the height of target spacecraft is assumed as 356 km in a circular orbit, and the mean angular velocity of the target spacecraft x = 0.0654°/s. The initial mass of the chaser is 200 kg and the initial mass of propellant is 50 kg. The size of thrusts are 100N in the x-axis direction, 90N in the y-axis direction and 80N in the z-axis direction, respectively. The massflow-rate of propellant of the chaser’s thrusters are 10 g/s in the x-axis direction, 9 g/s in the y-axis direction and 8 g/ s in the z-axis direction, respectively. According to Eq. (27), the rate of change of the line of nodes 2 X_ J 2  2:06071  1014 a7=2 ðcos iÞð1  e2 Þ :

ð95Þ

For a circular orbit of 356 km altitude and 52 deg inclination (International Space Station), the regression of nodes would be 4.789 deg/day. Solar radiation produces a force on a spacecraft in the Sun-spacecraft direction (Fehse, 2003): !

us : F SP ¼ qW~

ð96Þ

Per unit of mass of the satellite this is ~ cSP ¼ q

W ~ us ; m0

ð97Þ

where q is the radiation momentum flux, W is the cross section of the spacecraft, m0 is the initial mass of the chaser and ~ us is the Sun-spacecraft direction unity vector. The radiation momentum flux varies periodically with the orbit of the Earth around the Sun and is at aphelion : at perihelion :

q ¼ 4:38  106 N=m2 ; 6

2

q ¼ 4:68  10 N=m :

ð98Þ ð99Þ

Therefore gx, gy and gz are variables and they vary little during forming and keeping fast fly-around manoeuvres, then, in order to facilitate simulation, these coefficients can be seen as constants. Suppose that the sum of the perturbation and nonlinear factors caused by J2-effect and the solar radiation are gx = gy = gz = 0.00056N.

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431

ð100Þ where Tc is the time of the coplanar adjustment stage. Then according to Eqs. (27) and (28), the change of relative position and the relative velocity of the chaser in the x-axis during the coplanar adjustment are divided into the zerothrust stage and the decelerating stage as follows: In the zero-thrust stage zðtbz Þ ¼ 25:184 cos xtbz  61:1621 sin xtbz þ 74:816; vz ðtbz Þ ¼ 4 cos xtbz  1:647 sin xtbz :

1500

The relative coordinates of the three axes/m

Based on the vision measurement, the initial relative position and the initial relative velocity of the chaser are (1000m,  1000m, 100m) and (2m/s, 2m/s, 4m/s) and the target spacecraft is within the cone-shaped range of vision of the chaser. According to Eq. (23), the thrust acceleration in the z-axis is h i az ¼ 80 0:004 þ 1:28  107 T c þ 8:192  1012 ðT c Þ2 ;

1429

The change of x 1000

500

The change of z 0

−500

−1000

The change of y

−1500 0

50

100

150

200

250 time/s

300

350

400

450

500

Fig. 3. The change of x, y, z during the process of forming and keeping fast fly around.

ð101Þ ð102Þ

0.5

2

The change of ax,ay and az.(m /s)

0.4

In the decelerating stage 6

þ 2:38  10 tcz þ 66:21;

ð103Þ 6

0 ¼ 3:6 cos xtcz þ 2:4 sin xtcz þ 2:38  10 :

ð104Þ

The relative position and the relative velocity of the chaser are (937.5m, 937.5m, 0m) and (2m/s, 2m/s, 0) at the end of the coplanar adjustment stage. Based on the method presented in this paper, the time of the zero-thrust stage tbz = 18.75 s and the time of the decelerating stage tcz = 12.5 s during the coplanar adjustment. The time used to enter the fly-around orbit Te = 71.25s . In the process of the chaser enter the relative fast flyaround orbit the time of the zero-thrust stage tbx = 66.25 s and the time of the decelerating stage tcx = 5 s in the x-axis. The time of the zero-thrust stage tby = 63.8 s, the time of the decelerating stage tcy = 3.725 s and the time of the accelerating stage tay = 3.725 s in the y-axis. The relative position and the relative velocity of the chaser are (800 m,  800m, 0m) and (0m/s, 2m/s, 0m/s) when the chaser enter the fly-around orbit. The fast fly-around orbit is divided into 12 thrust arcs through the isochronous interpolation method and the time of each thrust arc is 30 s. There are three types of stages (the accelerating stages, the zero-thrust stages and the decelerating stages) in one thrust arc, and the magnitudes and the directions of thrust accelerations of the chaser in the x-axis and the y-axis are different in different thrust arc. Fig. 3 shows the change of x, y, z during the process of forming and keeping fast fly around. It can be seen from Fig. 3 that after the chaser enters the fast fly-around orbit, the relative position of the chaser in the z-axis keeps zero and the relative position of the chaser has periodical change in the x-axis and in the y-axis. Fig. 4 shows the change of ax and ay during the process of forming and keeping fast fly around. The magnitudes of

***: a

0.2

ooo: a

x

y

0.1 0

−0.1 −0.2 −0.3 −0.4 −0.5

0

50

100

150

200

250 time/s

300

350

400

450

500

Fig. 4. The change of ax, ay during the process of forming and keeping fast fly around.

8 6

ooo: Zero−thrust time ...: Accelerating time

4

The change of vx (m/s)

0 ¼ 36:816 cos xtcz  55:05 sin xtcz

0.3

***: Decelerating time 2 0 −2 −4 −6 −8

0

50

100

150

200

250 time/s

300

350

400

450

500

Fig. 5. The change of vx during the process of forming and keeping fast fly around.

the thrust accelerations in three axes are (0.4164 m/ s2, 0.3748 m/s2) at the end of the first fast fly-around cycle. Fig. 5 shows the change of vx and during the process of forming and keeping fast fly around. It can be seen from Fig. 5 that vx changes periodically in the keeping fast flyaround stage, and the time of each thrust arc will decreases

The change of vy and vz (m/s)

1430

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431 6

5. Conclusions

4

This paper presents a new switch control laws under constant thrust for the chaser fast flying around the target spacecraft along a specified fly-around trajectory. In this paper, the relative position parameters of the target spacecraft are obtained by using the vision measurement and the target fly-around positions are calculated through the isochronous interpolation method. The switch control laws are obtained based on the thrust acceleration sequences and the on time of the thrusters in three axes which can be computed by the time series analysis method. The perturbations and fuel consumptions are addressed during the computation of the on time of the thrusters. It is proved that, with the switch control laws, the chaser will fast fly around the target spacecraft along the specified fly-around trajectory.

2

0

−2

*** v z

−4

... vy

−6

0

50

100

150

200

250

300

350

400

450

500

time/s

Fig. 6. The change of vy, vz during the process of forming and keeping fast fly around.

with the increasing of the magnitude of the thrust acceleration in the x-axis. Due to the thrust acceleration is constant in the x-axis, therefore during the fast fly-around, the direction of the relative velocity in the x-axis has oscillatory behavior to keep the chaser fly around the target spacecraft. Fig. 6 shows the change of vy and vz during the process of forming and keeping fast fly around. It can be seen from Fig. 6 that vy changes periodically in the keeping fast flyaround stage, vz decreases to zero at end of the coplanar adjustment stage, then has no change (vz = 0) during the process of keeping fast fly around. Due to the thrust acceleration is constant in the y-axis, therefore during the fast fly-around, the direction of the relative velocity in the yaxis has oscillatory behavior to keep the chaser fly around the target spacecraft. Fig. 7 shows the change of the chaser’s trajectory during the process of forming and keeping fast fly around. It can be seen from Fig. 7 that there are 12 thrust arcs in the fast fly-around orbit. It is obvious that with the switch control laws the chaser can fly around the target spacecraft along the specified fly-around trajectory.

100 80

z

60 40

The trajectory of the chaser

20 0 1500 1000

1500

500

1000 0

500

0 −500 −500 −1000 −1000 y x −1500 −1500

Fig. 7. The chaser’s trajectory during forming and keeping fast fly around.

Acknowledgments This work is supported by the NSFC (60727002, 60774003, 60921001, 90916024), the MOE (20030006003), the COSTIND (A2120061303) and the National 973 Program (2005CB321902). References Enrico, C. Drag-free and attitude control for the GOCE satellite. Automatica 44, 1766–1780, 2008. Fehse, W. Automated rendezvous and docking of spacecraft. Cambridge Univ. Press, Cambridge, England, U.K.,, pp. 79–90, 2003. Geller, D.K. Linear covariance techniques for orbital rendezvous analysis and autonomous onboard mission planning. Journal of Guidance Control and Dynamics 29 (6), 1404–1414, 2006. Kukush, A., Huffel, S.V.. Consistency of elementwise weighted total least squares estimator in a multivariate errors-invariables model AX = B, vol. 59. Springer-Verlag, Metrika, pp. 75–97, 2004. Kojima, H. Fly-around motion control using geometrical differential theory, in: Tenth International Space Conference of Pacific-Basin Societies, Advances in the Astronautical Science, vol. 117, AAS-03383, 2003. Luo, J., Yang, Y., Yuan, J. Trajectory design and control of in-plane fast controlled flyaround. J. Astronaut. 11 (6), 1389–1392, 2006. Luo, J., Zhou, W., Yuan, J. A general method of trajectory design and guidance for fast satellite circumnavigation. J. Astronaut. 28 (3), 628– 632, 2007. Masutani, Y., Matsushita, M., Miyazaki, F. Fly around maneuvers on a satellite orbit by impulsive thrust control, in: Proceedings of the 2001 IEEE International Conference on Robotics Automation, Seoul Korea, pp. 421–426, 2001. Philip, N.K., Ananthasayanam, M.R. Relative position and attitude estimation and control schemes for the final phase of an autonomous docking mission of spacecraft. Acta Astron. 52, 511–522, 2003. Richards, A., Schouwenaars, T., How, J., et al. Spacecraft trajectory planning with avoidance constraints using mixed-integer linear programming. J. Guid. Control Dynam. 25 (4), 755–764, 2002. Schouwenaars, T., How, J.P., Feron, E. Decentralized cooperative trajectory planning of multiple aircraft with hard safety guarantees. AIAA Paper, pp. 2004–5141, 2004. Straight, S.D., Tragessor, S.G., Howard, R., et al. Maneuver Design for Fast Satellite Circumnavigation. Air Force Institute of Technology, 2004. Wu, H., Jiang, Y., Zhang, Z., et al. On-board propulsion technologies for microsatellites. J. Rock. Prop. 32 (3), 40–43, 2006.

Y. Qi, Y. Jia / Advances in Space Research 48 (2011) 1421–1431 Wang, F., Cao, X., Chen, X. On-orbit-servicing spacecraft flyaround orbit design and LQR keep control in eccentric orbits, in: Second IEEE Conference on Industrial Electronics and Applications, 2007. Wang, P., Yuan, J., Fan, J. Discussion on non-Keplerian orbit. J. Astronaut. 30 (1), 37–41, 2009. Wertz, J., Larson, W. Space Mission Analysis and Design. Microcosm and Kluwer Academic Publishers, Dordrecht, 1991. Yamanaka, K., Ankersen, F. New state transition matrix for relative motion on an arbitrary elliptical orbit. J. Guid. Control Dynam. 25 (1), 60–66, 2002. Yasuhiro, M., Motoshi,. M., Fumio, M. Flyaround maneuvers on a satellite orbit by impulsive thrust control, in: Proceedings of the 2001

1431

IEEE International Conference on Robotics Automation Seoul, Korea, May 21–26, 2001. Zachary, J.F. Orbit determination using modern filters/smoothers and continuous thrust modeling. Massachusetts Institute of Technology, 2008. Zhou, W., Yuan, J., Luo, J. Design and control of long-term flyaround trajectory for target spacecraft in non-coplanar elliptical orbit. Chin. Space Sci. Technol. 8 (4), 20–25, 2006. Zhu, Y. Trajectory Planning and Control for Spacecraft Proximity Relative Motion. National University of Defense Technology, Changsha, 2009.