Bounded Relative Orbits in the Zonal Gravitational Field for Formation Flying

Yanchao He,Ming Xu

DOI: https://doi.org/10.2514/1.g003107

2018-01-01

Journal of Guidance Control and Dynamics

Abstract:Open AccessEngineering NoteBounded Relative Orbits in the Zonal Gravitational Field for Formation FlyingYanchao He and Ming XuYanchao HeBeihang University, 100191 Beijing, People’s Republic of China*Ph.D. Student, School of Astronautics; .Search for more papers by this author and Ming XuBeihang University, 100191 Beijing, People’s Republic of China†Associate Professor, School of Astronautics; (Corresponding Author).Search for more papers by this authorPublished Online:3 Apr 2018https://doi.org/10.2514/1.G003107SectionsRead Now ToolsAdd to favoritesDownload citationTrack citations ShareShare onFacebookTwitterLinked InRedditEmail AboutNomenclatureasemi-major axisaz^normal acceleration of perturbations acting on the spacecraft in the local-vertical-local-horizontal (LVLH) frameeeccentricityh, hangular momentum vector and its magnitude (equal to |h|)i, i˙inclination and its derivative in terms of timeJlthe zonal harmonics coefficients of degree l (l≥2)Kkinetic energy of the spacecraftM, M˙mean anomaly and its derivative in terms of timePmap of the flow associated with relative motion equationsPl(sinφ)Legendre associated polynomialsRmean equatorial radius of the oblate bodyrposition vector of the spacecraftr, r˙distance of the spacecraft from the center of the body (equal to |r|) and the radial velocityUpotential of spacecraft subject to the attraction of an axisymmetric bodyu, u˙argument of latitude and its derivative in terms of timeμgravitational parameter of the oblate bodyφlatitude of the spacecraftρrelative position of the deputy with respect to the chief in the LVLH frameω, ω˙argument of perigee and its derivative in terms of timeλ, λ˙mean argument of latitude: ω+M, and its derivative in terms of timeΩ, Ω˙right ascension of the ascending node and its derivative in terms of timeω, ωx, ωy, ωzangular velocity of spacecraft in the LVLH frame and its componentsΔUdisturbing functionI. IntroductionIn recent years the studies on distributing the functionality of large spacecraft into a group of smaller and cooperating ones keep constantly emerging due to its successful application in formation-flying missions, such as TerraSAR-X (TSX) & TanDEM-X (TDX) [1], PRISMA [2], and GRACE [3]. In particular, the problem of searching naturally bounded relative orbits is of great significance for various distributed space missions, wherein the spacecraft remain within the relative position boundary while no or only little station-keeping effort is required [4,5].Much previous research has been dedicated to this problem from both analytical and numerical points of view since it was first proposed by Schaub and Alfriend [6]. They developed the first-order constraints for J2 bounded relative orbits, termed the standard J2-invariant bounded constraints, by linearizing the rate of differential mean right ascension of the ascending node (RAAN) and argument of latitude and setting them to be zero. The initial mean orbital elements setup of the relative orbits based on this condition can effectively minimize the relative drift, especially its along-track component over the certain periods of time. After this, studies around the bounded conditions and their improved versions have received broad attention from subsequent researchers. For example, Gurfil [7] obtained a more general J2-invariant orbits by using the nonosculating orbital elements under the concept of gauge-generalized equations. More recently, an additional correction factor was introduced to add into the standard invariant constraints enhancing the performance of J2 invariance in [8].Equivalently, the bounded relative motion can also be ensured by forcing the nodal period TΩ and drift of RAAN per nodal period ΔΩ of any two spacecraft to be equal on average. Noticing this fact, Xu et al. [9] presented a purely numerical method to search the J2-perturbed bounded relative orbits using the Poincaré surface of section and differential correction algorithm by improving Koon et al.’s work [10] on the investigation of the relative motion problem in the framework of dynamical systems theory. Similarly but furthermore, Baresi and Scheeres [11] studied naturally bounded relative orbits using quasi-periodic invariant tori calculated via stroboscopic maps, and they further extended their approach to explore the bounded relative orbits about celestial bodies such as Earth and the highly oblate primary of the binary asteroid when dealing with the zonal perturbations up to J5.To avoid the heavy numerical computations, several works about the analytical methods were also proposed by means of deriving the closed-form solutions of bounded relative motion. Two exact necessary conditions for boundedness of the relative motion were developed in [4], including the matching of the radial periods and the orbital angles, which was further used as a stepping stone for deriving single-impulse control of inter-spacecraft distance guaranteeing bounded relative motion [5]. Using the Cid’s intermediary and retaining the original structure of the full J2 Hamiltonian, Lara and Gurfil [12] found an integrable approximation for the relative motion problem to predict long-term bounded relative orbits with arbitrary inclinations. As the two quantities (TΩ and ΔΩ) are the essentials to bounded relative motion, Chu et al. [13] derived the analytical expressions of them to yield the bounded relative orbits with the action-angle variables instead of obtaining the complete solution of the J2-perturbed orbital motion.The multiple-shooting method, originally applied in the restricted three-body problem (RTBP) [14,15] and even the restricted four-body problem (RFBP) [16], has shown its success in finding the continuous quasi-periodic trajectories over long periods of time in the real ephemeris model. As the search of bounded relative orbits in our work is also the procedure of identifying the quasi-periodic orbit under higher-order perturbations with bounded constraints difficult to satisfy, it may also be beneficial to use this method to compute the initial value for bounded relative motion.Different from most of the previous works that model bounded relative motion under J2 perturbation, the present Note is devoted to the investigation of bounded relative orbits constraints considering the higher-order perturbations Jl (l>2) from the aspects of accurate modeling of relative dynamics around the Earth. First, to extend the standard J2-invariant bounded constraints, the proposed analytical method derives the first-order approximation of bounded constraints by initializing the chief on the frozen orbit under higher-order zonal perturbations, such as J4. This enables a better physical insight into the bounded relative motion compared with the traditional methods based only on J2 and the numerical approaches because, on the one hand, the proposed analytical one incorporates the long-periodic variations of the perturbations that the traditional methods neglected; on the other hand, it deals with the generation of bounded relative motion based on the orbital elements and the bounded constraints, in accordance with habit of the orbit design engineers who prefer the classical orbital elements than other quantities such as hybrid variables (r,r˙,h,i,u)T employed later in this Note. Then, to remedy the shortcoming of the analytical approach, the bounded relative orbits are yielded via numerical procedure that computes the quasi-periodic orbit by mean of the multiple-shooting approach. As the methods presented in this work directly exploit the higher-order perturbations environments, the accurate modeling enables a wider application of them in the formation-flying mission about the oblate body, compared with the standard J2 constraints and their improved forms.II. Equations of Relative Motion Under Higher-Order Zonal Harmonic PerturbationsTo analyze the bounded relative orbits of formation flying, it is necessary to develop the exact dynamics of the relative motion in an axially symmetric gravitational field of an oblate body rather than adopt the usual and simpler dynamics of only second order or its linearized form. Thus, the gravitational potential considering the higher-order zonal harmonics is given by [17] U=μr+ΔU=μr[1−∑l=2∞Jl(Rr)lPl(sinφ)](1)where μ is the gravitational parameter, R is the mean equatorial radius of the oblate body, r is the distance of the spacecraft from the center of the body, φ is the latitude, Jl are the zonal harmonics coefficients, and Pl(sinφ) are the Legendre associated polynomials of degree l. In the following parts, all the length and time values are scaled by the characteristic length R and time R3/μ, respectively. Therefore, both the equatorial radius and gravitational parameter are equal to 1.In the potential determination process, the order of potential depends on what degree of Legendre polynomials is expanded up to. The derivation procedure with zonal terms of nonspherical gravitational disturbing function up to the arbitrary Jl (l>2) order maybe a bit tedious and trivial, but the detailed expressions with terms from J2 up to J15 have already been addressed in [18].For the purpose of describing the exact relative motion of spacecraft attracted by the nonspherical gravity of an oblate body, the body-centered initial (BCI) coordinate frame and the chief-fixed local-vertical-local-horizontal (LVLH) rotating coordinate frame are defined, respectively, as shown in Fig. 1. First, the origin of BCI frame is located at the central body with the axis X^ pointing to the vernal equinox, axis Z^ along the rotation axis of the body, and axis Y^ completing the right-handed frame. Second, as for the LVLH frame, its origin coincides with the mass center of the chief, the axis x^ is directed from the chief radially outward, the axis z^ is along the direction of instantaneous angular momentum of the chief, and the axis y^ is finally orthogonal to the previous ones also satisfying the right-handed rule. It needs to be pointed out that, unless particularly noted, any variable associated with subscript “c” indicates that it belongs to the chief, whereas variable with “d” indicates that of the deputy.Fig. 1 BCI and LVLH coordinate frames.Let ic be the inclination, Ωc the RAAN, and uc the argument of latitude; rc refers to the distance of the chief in the BCI frame, r˙c refers to the radial velocity of rc, and hc refers to the angular momentum. Considering the zonal harmonics of gravity acting on the spacecraft, the state of the chief can be described by the variables (rc,r˙c,hc,ic,uc)T [19]. Thus, they can be written by the following equations as r¨c=−1r2+h2r3+∑l=2∞(l+1)Jlrc−(l+2)Pl(sinicsinuc)(2a)h˙c=−∑l=2∞Jlrc−(l+1)∂Pl(sinicsinuc)∂uc(2b)i˙c=rccosuchcaz^(2c)u˙c=hcrc2−Ω˙ccosic(2d)Ω˙c=rcsinuchcsinicaz^(2e)az^=−1rcsinuc∑l=2∞Jlrc−(l+1)∂Pl(sinicsinuc)∂ic(2f)where az^ is the normal acceleration of perturbations acting on the chief and expressed in the LVLH frame.The angular velocity of chief in the LVLH frame is given by ωc=Rz(uc)Rx(ic)[00Ω˙c]+Rz(uc)[i˙c00]+[00u˙c](3)where Rx(ic) and Rz(uc) are the elementary rotation matrixes around the axes X^ and Z^ by the corresponding angles, respectively. Thus, the components of the angular velocity are ωcx=i˙ccosuc+Ω˙csinucsinic(4a)ωcy=Ω˙ccosucsinic−i˙csinuc(4b)ωcz=u˙c+Ω˙ccosic(4c)Substituting Eq. (2) into Eq. (4) yields the following equations: ωcx=−1hcsinuc∑l=2∞Jlrc−(l+1)∂Pl(sinicsinuc)∂ic(5a)ωcy=0(5b)ωcz=u˙c+Ω˙ccosic=hcrc2(5c)Therefore, the components of angular acceleration are obtained by taking derivative of Eq. (5) in terms of time [20]: ω˙cx=−1hc2sinuc(∑l=2∞Jlrc−(l+1)∂Pl(sinicsinuc)∂ic)⁢(∑l=2∞Jlrc−(l+1)∂Pl(sinicsinuc)∂uc)+r˙chcsinuc(∑l=2∞(l+1)Jlrc−(l+2)∂Pl(sinicsinuc)∂ic)+∂ωcx∂uc(hcrc2+cosichcsinic∑l=2∞Jlrc−(l+1)∂Pl(sinicsinuc)∂ic)−∂ωcx∂iccosuchcsinuc∑l=2∞Jlrc−(l+1)∂Pl(sinicsinuc)∂ic(6a)ω˙cy=0(6b)ω˙cz=ddt(hcrc2)=h˙crc2−2hcr˙crc3=h˙crc2−2ωczr˙crc=−∑l=2∞Jlrc−(l+3)∂Pl(sinicsinuc)∂uc−2hcr˙crc3(6c)Consequently, given the state of the chief (rc,r˙c,hc,ic,uc)T, the angular velocity ωc and its corresponding angular acceleration ω˙c can be obtained via Eqs. (5) and (6).Consider the relative motion of formation flying in the LVLH frame. The position and velocity of the chief (deputy) are denoted by rc (rd) and vc (vd), respectively. The relative position of the deputy with respect to the chief in the LVLH frame is represented by ρ=rd−rc=(x,y,z)T. Hence, the relations between the above variables can be expressed as rd=rc+ρ(7)vd=vc+ω˙×rd+ρ˙(8)The position of the deputy expressed in the BCI frame is given by rId=Rz(−Ωc)Rx(−ic)Rz(−uc)(rc+ρ)(9)Thus, the distance from the center of the body to the deputy is obtained as rd=(rc+x)2+y2+z2(10)and the third component of rId expressed in the BCI frame is rdz=(r+x)sinicsinuc+ysiniccosuc+zcosic(11)Then, the gravitational potential of the deputy is given by Ud=−1rd+∑l=2∞Jlrd−(l+1)Pl(rdzrd)(12)and the kinetic energy is Kd=12vdTvd=12(x˙+r˙−yωcz)2+12[y˙+(r+x)ωcz−zωcx]2+12(z˙+yωcx)2(13)Based on the Lagrangian formulation [19] of relative motion of the deputy in the LVLH frame shown as follows: ddt(∂Kd∂ρ˙)−∂Kd∂ρ+∂Ud∂ρ=0(14)Substituting Eqs. (12) and (13) into Eq. (14) and making some manipulations yields the following nonlinear equations of relative motion in higher-order zonal harmonic perturbations x¨=2y˙ωcz+xωcz2+yω˙cz−zωcxωcz+rωcz2−r¨−∂Ud∂x(15a)y¨=−2x˙ωcz+2z˙ωcx−xω˙cz+y(ωcx2+ωcz2)+zω˙cx−2r˙ωcz−rω˙cz−∂Ud∂y(15b)z¨=−2y˙ωcx−xωcxωcz−yω˙cx+zωcx2−rωcxωcz−∂Ud∂z(15c)To conclude, the state of the chief (rc,r˙c,hc,ic,uc)T can be determined by Eq. (2). After this, the angular velocity ωc=(ωcx,ωcy,ωcz)T and its corresponding angular acceleration ω˙c=(ω˙cx,ω˙cy,ω˙cz)T are naturally computed via Eqs. (5) and (6). Then, according to Eqs. (10–12), the gradient of gravitational potential of the deputy ∂Ud/∂ρ=(∂Ud/∂x,∂Ud/∂y,∂Ud/∂z)T can be calculated; likewise, the gradient of kinetic energy ∂Kd/∂ρ=(∂Kd/∂x,∂Kd/∂y,∂Kd/∂z)T can be obtained from Eq. (13). Having known all the quantities mentioned above, the relative motion of the formation flying is presented according to Eq. (15), which lays the basis of Sec. IV.III. Bounded Relative Orbits ConstraintsAccurate formulation of initial conditions to set up the relative orbits is essential to the orbit design of formation flying. To this end, one should take into account the higher-order perturbations to develop the conditions for the generation of bounded relative orbits. Assume a scenario of two-spacecraft formation in which one is assigned to be the chief and the other is the deputy. The orbital elements (including semi-major axis, eccentricity, inclination, argument of perigee, RAAN, mean anomaly, and mean argument of latitude) of the chief (deputy) are correspondingly denoted by ac (ad), ec (ed), ic (id), ωc (ωd), Ωc (Ωd), Mc (Md), and λc=ωc+Mc (λd=ωd+Md), and δ(⋅)=(⋅)d−(⋅)c are referred to as the difference between the chief and deputy. To remark the distinction with Sec. II, it should be noted here that all the orbital elements in this section indicate the mean ones, which will be used to develop the bounded relative orbits constraints under the averaged perturbations over the mean anomaly M from 0 to 2π, that is, ΔU¯=(1/2π)∫02πΔU dM.Following the idea of J2-invariant orbit proposed by Schaub and Alfriend [6], the drift rates of the mean argument of latitude and RAAN are the two key factors for the design of bounded relative orbits; it is required that the differences of them (δλ˙=δΩ˙=0) be zero. The initial setup of orbital elements for the formation satisfying this requirement can maintain the two spacecraft on orbits synchronously with no along-track drift and thus be bounded.However, when higher-order (Jl, l>2) terms of perturbations are taken into consideration, the case for bounded relative orbits will be much different from that under only J2 perturbation. According to the Lagrange planetary equations [21] substituted into the averaged perturbing function, if merely the J2 perturbation of gravitational potential is included, the semi-major axis, eccentricity, and inclination keep constant on average, whereas if higher-order zonal perturbations are considered [18], the long-periodic value ω begins to take effect in the perturbations and only the value of mean semi-major axis is a constant. As expressed by the Lagrange planetary equations, the secular variations of eccentricity e˙, inclination i˙, argument of perigee ω˙, RAAN Ω˙ and mean anomaly λ˙ depend on semi-major axis (a), eccentricity (e), inclination (i), and argument of perigee (ω). Thus, it has to be mentioned that the variations of these orbital elements (e, i, and ω) are the obstacles to design the bounded relative orbits. Thus, it needs to remove the variation of some elements. Fortunately, e˙=i˙=0 can be naturally satisfied by setting ω to be either 90° or 270°, which requires that the orbit should be on the frozen state. Thus, with a required semi-major axis and the eccentricity given, the corresponding inclination can be computed by using any standard numerical method according to ω˙=0. Hence, to accommodate the higher-order perturbations effect, it is necessary to imposing the following constraints to guarantee the bounded relative motion: δe˙=∂e˙c∂acδa+∂e˙c∂ecδe+∂e˙c∂icδi+∂e˙c∂ωcδω=0(16a)δi˙=∂i˙c∂acδa+∂i˙c∂ecδe+∂i˙c∂icδi+∂i˙c∂ωcδω=0(16b)δλ˙=∂λ˙c∂acδa+∂λ˙c∂ecδe+∂λ˙c∂icδi+∂λ˙c∂ωcδω=0(16c)δΩ˙=∂Ω˙c∂acδa+∂Ω˙c∂ecδe+∂Ω˙c∂icδi+∂Ω˙c∂ωcδω=0(16d)Based on the Lagrange planetary equation, as the chief is on the frozen orbit, e˙c=i˙c=0, Eqs. (16a) and (16b) will deliver that δω=0, and so the bounded constraints are reduced to δλ˙=∂λ˙c∂acδa+∂λ˙c∂ecδe+∂λ˙c∂icδi=0(17a)δΩ˙=∂Ω˙c∂acδa+∂Ω˙c∂ecδe+∂Ω˙c∂icδi=0(17b)Therefore, given the deviation of inclination δi, Eq. (17) can compute the differences in mean semi-major axis, and mean eccentricity as follows: δa=A2B3−A3B2A1B2−A2B1δi(18a)δe=A3B1−A1B3A1B2−A2B1δi(18b)where Ai, Bi (i=1,2,3) denote auxiliary variables, defined by A1=∂λ˙c∂ac,A2=∂λ˙c∂ec,A3=∂λ˙c∂ic(19a)B1=∂Ω˙c∂ac,B2=∂Ω˙c∂ec,B3=∂Ω˙c∂ic(19b)And the rate of mean orbital elements (e˙c, i˙c, λ˙c, and Ω˙c) as well as the auxiliary variables (A1–A3, B1–B3) are derived from Lagrange planetary equations with the averaged zonal perturbations (refer to Appendix). Consequently, given the prescribed inclination difference between the deputy and the chief δi, the differences in semi-major axis and eccentricity can be computed from Eq. (18).By setting the initial conditions the relative orbits with the osculating orbital elements based on the mean orbital elements acquired above after the mean-to-osculating transformation [21,22], we can achieve the bounded relative motion with the minimum secular drift caused by the effects from higher-order perturbations.As for the Earth case, the first three main gravitational harmonics (corresponding to the zonal harmonic coefficients J2, J3, and J4) of the geopotential are taken into consideration. Of course, in our method the order of zonal harmonic perturbations can be expanded up to higher than J4. In this case, the initial conditions assigned to the chief and the deputy are determined by Eq. (18), as displayed in Table 1.As discussed earlier, to ensure the bounded relative motion, the chief should be on the frozen orbit under zonal harmonic perturbations. Figure 2 exhibits the variations of the osculating (solid line) and mean (dashed line) orbital elements in the time window of interest, which validates the frozen property of the chief. This is the basic guarantee to make the bounded relative motion possible in our work.Fig. 2 Osculating and mean orbital elements.With the setting of initial relative motion geometry in Table 1, Fig. 3 shows the evolution of J4 bounded relative orbits in the chief’s LVLH frame over 5 days (equivalent to about 73 orbital revolutions); correspondingly, the time histories of the relative distance and along-track distance are plotted in Fig. 4.Fig. 3 J4 bounded relative orbits in the LVLH coordinate frame.Fig. 4 a) Relative distance and b) along-track distance varying with time.As can be seen from Figs. 3 and 4, the relative configuration of chief and deputy is maintained as bounded as the drifts in relative distance and along-track motion can be effectively suppressed. These show that the initial conditions computed with the proposed analytical method are valid for yielding the naturally bounded relative orbits under the higher-order zonal harmonic perturbations.However, because the initial conditions for bounded relative orbits are derived analytically based on the first-order linearization approximation to the bounded constraints, it is necessary to verify the accuracy and reasonability of this approach. Thus, the validity of the linearized equations is performed by the comparison of relevant quantities (δe˙, δi˙, δλ˙, and δΩ˙) involved in the derivation of bounded constraints for the actual values and their approximations, as shown in Fig. 5. It is easy to observe from the plots that the linearized approximations of relative eccentricity rate and the relative inclination rate are in good agreements with the actual ones, whereas the linearizations to approximate the actual values of the relative RAAN rate and the relative mean argument of latitude are with certain differences below the order of magnitude 10−12, which is anyway a tiny value. Therefore the linear approximation is valid for the motion of small size of relative orbits in the short time. Of course, like the previously devised J2-invariant relative orbits [6], it is due to the unavoidable linearization error that the proposed analytical method would be destined to be inapplicable in the long-term and large-size relative motion case.Fig. 5 The comparison of the actual rates of relative orbital elements with their linear approximations.To show the improvements of our bounded conditions compared with the J2 condition, we initialize the chief on the same orbit, but the deputy’s orbits are obtained from ours and J2 method independently. The behaviors of relative distances starting from the J4 (ρJ4) and J2 (ρJ2) bounded conditions, as well as their difference (ρJ4−ρJ2), are provided by Fig. 6. As manifested in Fig. 6c, the difference has a downtrend tendency, around 3.5 m in 20 days, which shows that the relative orbit from the J2 bounded conditions drifts quicker than that from the J4 conditions under the same higher-order perturbations. This tells that our method can guarantee the better precision to achieve the bounded relative motion compared with the traditional J2 one.Fig. 6 The relative distance difference between J4 and J2 bounded conditions.In addition, as the bounded conditions are studied under the zonal harmonic perturbations, to make them more applicable in the practical perturbations environments, it is interesting to investigate the effect of the nonzonal gravitational perturbations, such as the tesseral harmonic terms (J22) on the bounded relative orbits. Figure 7 shows the drift of relative distance due to the addition of J22 term relative to the relative distance only in the zonal perturbations. With the time elapsed, the drift is magnified gradually, which can arrive at around 0.12 km in 20 days.Fig. 7 The relative distance drift caused by tesseral harmonic terms (J22).It has to be mentioned that even though the bounded relative orbits have been achieved by the proposed method, the bounded configuration will not be held any more if the expected time interval is extended to a much longer duration, such as 30 days. As illustrated by Fig. 8, the relative orbit initialized as the same orbital elements setup with Table 1 forms an closed distorted band, which indicates that the deputy drifts away from the chief; correspondingly the relative distance is no longer bounded with a drift approximately 0.36 km during this period.Fig. 8 Relative orbits in the LVLH frame in one month: a) the three-dimensional relative motion; b) the along-track distance; and c) relative distance.As pointed out by some previous researchers [12,23], because of the errors lying in the linearization to the bounded constraints and the mean-to-osculating transformation [21,22] on the orbital elements, the proposed analytical approach only succeeds to hold the stable bounded relative geometry in a very short time and small size scale. Therefore, we need to switch to a numerical approach for the investigation of the bounded relative motion.IV. Bounded Relative Quasi-Periodic Orbits by the Multiple-Shooting ApproachBecause of the drawback of the analytical method to find the bounded relative orbits with slight along-track drift, this section is dedicated to the generation of quasi-periodic orbit that satisfies the bounded constraints in a fully numerical way. To this end, our attention is focused on how to achieve the bounded relative motion based on the selected initial guess. The search of bounded relative orbits is actually a process of finding the quasi-periodic orbits, which meet the bounded constraints, to guarantee the bounded relative motion for the formation-flying spacecraft. Motivated by the multiple-shooting approach for identifying quasi-periodic solutions in real ephemeris model [14,15], in order to obtain the bounded relative orbits, the initial guess of the approach would be some kind of orbit as long as it can guarantee the bounded relative motion. It is worth noting that the initial orbit is not necessary to be the real one that happens in the perturbed physical environments, but should be strictly kept bounded in relative distance, which is because the resulting orbit by the multiple-shooting approach is very close to the initial input and thus highly depended on its boundedness. Basically, the relative periodic orbits can be the good options for the initial input of the numerical algorithm. Depending on how eccentric the reference orbit of the chief is, the Hill-Clohessy-Wiltshire (HCW) [21] and/or Tschauner-Hempel (TH) [24,25] solutions are qualified enough for the algorithm to iteratively compute the bounded relative quasi-periodic orbits, very close to the initial guess, under the complex perturbations.The whole time interval from t1 to tN is split into a series of evenly spaced subinterval with the number of (N−1), thus we can let the initial time point be t1, the final time point be tN, and the length of the subinterval Δt=tk+1−tk, k=1,2,…,N−1. The initial conditions to be input into the iterative procedure are Q=(x1,y1,z1,x˙1,y˙1,z˙1,…,xk,yk,zk,x˙k,y˙k,z˙k,…,xN,yN,zN,x˙N,y˙N,z˙N)T=(Q1,…,Qk,…,QN)Twith the total number of variable components 6N. Define a map P by the flow associated with relative motion equations [Eq. (15)] of formation flying after the time interval Δt. Let P(Qk) be the image of the point Qk. As all the points are on the same orbit, to find the value of Q, the condition can be imposed by Fk(Q)=P(Qk)−Qk+1=0 in which Fk denotes the kth block of six scalar equations. If we write all the conditions in a concise vector form, they can be expressed by F[Q1Q2⋮QN]=P[Q1Q2⋮QN−1]−[Q2Q2⋮QN]=0(20)where F denotes the vector of 6(N−1) equations to be solved by the variable Q with 6N components. Given a reliable initial value of variable Q(0), the jth iterate for the improved value Q(j) can be obtained via the Newton’s method until the convergence; that is, ‖F(Q(j−1))‖ is below the given threshold, such as 10−12: DF(Q(j−1))(Q(j)−Q(j−1))=−F(Q(j−1))(21)in which the differential matrix DF is formed by the differential of the map and the identity matrix as DF=[A1−IA2−I⋱⋱AN−1−I](22)In Eq. (22), Ak, k=1,2,…,N−1, stand for the 6×6 matrix DP(Qk), that is, the differential of the map evaluated at the point Qk during each iterate procedure. It is worth noting that at each step the equations number 6(N−1) is less than the unknown components of Q6N. Here, an additional condition that the Euclidean norm of consecutive iterative value ‖ΔQ(j)‖=‖Q(j)−Q(j−1)‖ is required to be minimum should be attached to the iterative procedure, which can be obtained by introducing the Lagrange function: L(ΔQ,λ)=ΔQTΔQ+λT(F(Q)+DF(Q)ΔQ)(23)where λ is the multiplier vector. And we can find the optimum solution ΔQ(j)=−DF(Q(j−1))T[DF(Q(j−1))DF(Q(j−1))T]−1F(Q(j−1))(24)Thanks to the methods of matrix decomposition and recursions proposed by [14], we can get rid of the huge burden of computing the inverse of sparse matrix [DF(Q(j−1))DF(Q(j−1))T]−1 directly. If we define a symmetric matrix M=DF(Q(j−1))DF(Q(j−1))T=[I+A1A1T−A2T−A2I+A2A2T−A3T⋱⋱⋱−AN−2I+AN−2AN−2T−AN−1T−AN−1I+AN−1AN−1T](25)and introduce the new variables Z(j−1) by Z(j−1)=M−1F(Q(j−1)), Eq. (24) is accordingly replaced by ΔQ(j−1)=−DF(Q(j−1))TZ(j−1)(26)To derive the expressions of Z and ΔQ, the following factorization on the symmetric matrix M is performed as M=[IL2I⋱⋱LN−1I][D1D2⋱DN−1]⁢[IL2TI⋱⋱LN−1TI](27)where the corresponding relations are yielded as {D1=I+A1A1TLk=−AkDk−1−1,k=2,3,…,N−1Dk=I+AkAkT−LkDk−1LkT,k=2,3,…,N−1(28)With the help of the relations above, the variable Z can be thus obtained from the following recursions: {X1=F1(Q(j−1))Xk=Fk(Q(j−1))−LkXk−1,k=2,3,…,N−1YN−1=DN−1−1XN−1Yk=Dk−1Xk−Lk+1TYk+1,k=N−2,N−3,…,1(29)in which X and Y are the auxiliary variables, with block components X1,…,XN−1, Y1,…,YN−1, and Y=Z(j−1). Consequently, the value of ΔQ(j) for iteration is acquired by Eq. (26).The procedure of performing multiple-shooting method for finding the quasi-periodic orbits is illustrated in Fig. 9, in which a) describes the correction between the neighboring nodes and b) is a global visualization of this algorithm. After the algorithm converges, in each subinterval, the resulting orbit is the closest agreement with the initial guess. To be specific, the resulting quasi-periodic orbit is depended on the initial guess, which can be anything as long as it guarantees boundedness in relative distance; furthermore, as for the resulting orbits in these subintervals, the terminal point of one orbit and the starting point of t