Relative Equilibrium and Stability of a Three-body Tethered Satellite in an Elliptical Orbit

This work focuses primarily on the relative equilibrium and stability of a three-body tethered satellite which center of mass moves along an elliptical orbit with arbitrary eccentricity in the central Newtonian field. We only consider the case with straight tethers, and the lengths of them are small compared to the size of the orbit. It means that the attitude motion of the satellite has no impact on its orbital motion. Additionally, it assumes that the lengths of the tethers are controllable. The law of the tether length control is determined on the basis of the dynamical model of the system so that there exists an Earth-pointing relative equilibrium of the tethered satellite. By analyzing the Floquet multipliers of the linearized system, Lyapunov instability regions and linear stability regions of the relative equilibrium are obtained in a plane of dimensionless parameters. Further nonlinear stability analysis is presented in the linear stability regions by means of a normal form technique and stability theorems on the Hamiltonian system. In particular, stability of nonresonant cases as well as fourth-order resonant cases are investigated in detail.


Introduction
Over the past few decades, the dynamics and control of tethered satellite systems have attracted a significant amount of attention. A tethered satellite system is comprised of a mother satellite and several subsatellites, which are connected with each other by thin and long tethers. There are plenty of important application prospects in the space solar power station, deep space exploration, space debris active removal and many other space missions [1,2,3,4,5,6].
In the investigation of attitude dynamics and control of tethered satellite systems, the satellite is regarded as a particle and the tether is simplified as a massless rigid rod. The one-piece dumbbell model and two-piece dumbbell model are the simplest and most commonly used simplified models. In the early days, Matteis [7] used the one-piece dumbbell model to study the effects of orbital eccentricity on the stability of a tethered satellite system, and it was pointed out that the stability criterion of elliptical orbit with small eccentricity was the same as that of circular orbit. Cui et al. [8] studied the stability for the deployment and retrieval of a tethered satellite system in an elliptical orbit by using the one-piece dumbbell model and the conditions for the existence of an equilibrium state were also obtained. Several years later, Celletti and Sidorenko [9] studied the attitude dynamics of the dumbbell satellite in the central gravitational field and noted that when the center of mass moves on a circular orbit, one can find a stable relative equilibrium, i.e. the so called local vertical. In the case of an elliptical orbit, however, there is no stable equilibrium even for small values of eccentricity. Additionally, Modi and Brereton [10,11] investigated periodic solutions and stability of a gravitygradient-oriented satellite using numerical method and pointed out that in an elliptical orbit which eccentricity is smaller than 0.35, the stable libration of a tethered satellite is possible.
In an elliptical orbit, even the equations of motion of the dumbbell model of a tethered satellite system become strongly nonlinear and then complex nonlinear phenomena appear. In this case, due to changes of the gravity gradient and orbital angular velocity, there is no equilibrium state and the tethered system begins tumbling. Using the one-piece dumbbell model, Hironor et al. [12] investigated the nonlinear dynamical behavior of a tethered satellite which center of mass moves in a circular orbit and an elliptical orbit, respectively. Misra et al. [13] studied the chaotic motion of a two-body tethered satellite, which is simplified as a dumbbell model, using numerical tools such as phase portraits, spectral analysis, Poincaré sections and Lyapunov exponents. Sideorenko and Celletti [14] studied the stability and bifurcation of periodic motions of a spring-mass model of tethered satellite subjected to both in-plane and outof-plane perturbations. They noted that the stable periodic motions only take place in the case of orbital eccentricity is sufficiently small. Jung et al. analyzed the nonlinear dynamical behavior of a tethered satellite system with a moving mass [15] and a three-body tethered satellite system undergoing deployment and retrieval [16] by using the two-piece dumbbell model.
Attitude control and vibration suppression of tethered satellites in an elliptical orbit have been an attractive topic due to the complex dynamic characteristics. The tether length control is one of the most popular control strategies of tethered satellite systems. Takeichi et al. [17] presented fundamental strategies for libration control and deployment of one-piece dumbbell model of a tethered system in an elliptical orbit, in which the control force was provided by a thruster carried by the subsatellite. Kojima et al. [18] studied the control method of librational motion of the two-piece dumbbell model of tethered satellite in an elliptical orbit, in which the delayed feedback control method using thrusters on subsatellites was used to stabilize the chaotic, librational motion of the satellite to a periodic motion. Kumar et al. [19] designed an open-loop control law to regulate the tether length to control the libration of a tethered system in an elliptical orbit. Williams [20] designed a feedback controller and adopted the method for adjusting the length of the tether to control the librational motion of a tethered satellite in an elliptical orbit into periodic motion and stability of the periodic solutions was analyzed using Floquet theory. Recently, Shi et al. [21] simplified a tethered satellite with a moving climber in a circular orbit into the two-piece dumbbell model and pointed out that a control strat-egy based on the sliding mode control is effective in suppressing the libration of the tethered system with tension control only. Yu et al. [22] used the dumbbell model to study the chaotic behavior of a tethered satellite induced by atmospheric drag and the oblateness of the earth and the chaotic motion was suppressed by a tether length control based on a sliding mode controller. Furthermore, Kojima et al. [23] developed an optimal acceleration climber transit control scheme for the two-piece dumbbell model of a tethered satellite in a circular orbit and the optimal control law was verified by experiments. Shi et al. [24] presented a control strategy to suppress the libration of the two-piece dumbbell model of a partial space elevator by the control torque of the climber generated by reaction wheels.
Investigations of stability of the tethered satellite system is necessary, since stability affects the success of the space mission. Kojima and Sugimoto [25] analyzed the stability of in-plane and out-of-plane motion of electrodynamic tethered systems in elliptic and inclined orbits. Stability analysis based on Floquet theory confirmed that a delayed feedback control with appropriate control gains could stabilize the librational inplane and out-of-plane periodic solutions. He et al. [26] studied the stability of a dumbbell model of a tethered satellite system in an elliptical orbit. Based on the linear system stability criteria, it was noted that to stabilize the tethered satellite system during deploying, station-keeping and retrieving, the eccentricity must satisfy e ∈ [0, 0.75), and some extra control force should be acted on the system. Using the one-piece dumbbell model, Zhang et al. [27] presented several criteria on the existence of periodic solutions for a tethered satellite system in an elliptical orbit and the uniqueness of periodic solutions in the circular orbits. Furthermore, the Lyapunov stability method and Barbashin-Krasovski theorem were used to analyze the stability of the equilibrium states of the tethered satellite system.
The gravity gradient can be used to stabilize a satellite in circular orbit. There is, however, no relative equilibrium in an elliptical orbit. As a consequence, an active attitude control should be considered to stabilize the satellite to Earth-pointing equilibrium. A method based on the regulation of inertia tensor to achieve attitude control in an elliptical orbit has been paid much attention. Connell [28] proposed a control strategy for an Earth-pointing satellite in an elliptical orbit by telescoping two identical masses outward from the satellite along the yaw axis. Additionally, Floquet theory was used to analyze the stability of relative equilibrium. Aslanov and Bezglasnyi [29] studied the plane motion of an axisymmetric satellite with a movable mass and the stability of two opposite relative equilibrium were investigated. More than that, using the swing-by technique, a continuous control law of the movable mass was constructed to stabilize the axis of symmetry of the satellite to the local vertical and to reorient the satellite from one stable equilibrium to the other. Markeev [30] studied a rigid satellite carrying a movable mass point in a central Newtonian gravitational field. The law of the relative motion of the mass point was obtained, which can maintain the relative equilibrium of the satellite in an elliptical orbit, as well as the stability of relative equilibrium was examined. Kholostova [31] studied the same problem and a nonlinear stability analysis of the relative equilibrium was performed based on the normalization method of area-preserving mapping and Hamiltonian system stability theory. Recently, based on the Hamiltonian approach, Burvo et al. [32] studied the dynamics of a rigid spacecraft with variable mass distribution. Several relative equilibrium and the corresponding control laws were derived in an elliptical orbit and the stability of relative equilibrium was investigated in detail.
In this paper, a three-body tethered satellite system in an elliptical orbit is simplified as the two-piece dumbbell model. To obtain an Earth-pointing relative equilibrium in an elliptical orbit with arbitrary eccentricity, the mass distribution of the system is controlled by regulating the length of the tether. The current work studies the control law of the tether length which makes the satellite existing a relative equilibrium in an elliptical orbit with arbitrary eccentricity and investigate the stability of relative equilibrium in detail. Previous studies on stability of tethered satellites are mainly based on Floquet theory, which estimates the stability by calculating the Floquet multiplier of the linearized system. We focus on additional nonlinear stability by using a method based on normal form technique and KAM theory [33,34,35,36], and stability of the nonresonant case and fourth-order resonant case is discussed in detail.
The rest of this paper is structured as follows. In Sect. 2, based on the two-piece dumbbell model, equations of motion of a three-body tethered satellite in an elliptical orbit with arbitrary eccentricity are derived. In Sect. 3, the control law of tether length is obtained, which ensures the relative equilibrium of the satellite in an elliptical orbit. In Sect. 4, the Hamiltonian of controlled system subjected to the planar perturbation is derived. In Sect. 5, coefficients of the normal form of Hamiltonian are calculated by means of an algorithm of constructing symplectic mapping. In Sect. 6, stability of the relative equilibrium of the three-body tethered satellite is examined, including the cases: mother satellite is in the lowest orbit, intermediate orbit and highest orbit. Finally, Sect. 7 is the summary and conclusion of this work.

Equations of motion
Consider a three-body tethered satellite system which is simplified as the classical two-piece dumbbell satellite, as shown in Fig. 1, in which O e is the center of attraction and O c is the center of mass of the system. Assume that the center of mass of the system moves along an elliptical orbit with arbitrary eccentricity. Let S c = {i 0 , j 0 } be the orbital coordinate frame which origin coincides with the center of mass O c , the axis i 0 is directed along the position vector r c of the center of mass with respect to the center of attraction O e to the center of mass O c , the axis j 0 is directed along a transversal to the orbit of the center of mass O c . In addition, we introduce another two local coordinate frames S 1 = {i 1 , j 1 } and S 2 = {i 2 , j 2 }, the axis i 1 is directed along the vector from mass point m 1 to m 2 , the axis i 2 is directed along the vector from mass point m 2 to m 3 . The angle between the unit vector i 0 and i 1 is expressed as ϕ 1 while ϕ 2 represents the angle between the unit vector i 0 and i 2 . The mass points m 1 , m 2 and m 3 represent the mother satellite and subsatellites. The position vector of mass points m 1 , m 2 and m 3 with respect to center of mass O c are represented as r 1 , r 2 and r 3 , respectively.
Furthermore, we make the following assumptions: (i) The mother satellite and subsatllites are point masses and the mass of tether is ignored; (ii) Only motions with straight strained tethers are considered, therefore, the tether is regarded as a massless rigid rod; (iii) The geometry of the satellite is much smaller than the radius of orbit hence the attitude motion of the satellite does not affect its orbital motion; (iv) Only the gravity gradient of the center of attraction is considered, and other perturbation torques, such as gravitational force of other celestial bodies, atmospheric drag and solar pressure, are ignored; (v) The tethers length L 1 and L 2 are controllable which are twice continuously differentiable functions of time with initial position ϕ 1 = ϕ 2 = 0.
According to the assumptions (ii) and (iv), the center of mass of the system is fixed at the origin of the orbital coordinate frame, i.e.
In addition, as shown in Fig. 1, there are the following relationships between the position vectors of the parti-cles m 1 , m 2 and m 3 and Therefore, the position vectors of particles m 1 , m 2 and m 3 with respect to the center of mass O c can be represented as in which The position vector of particle m j (j = 1, 2, 3) with respect to the center of attraction O e can be expressed as follows in which the position vector of the center of mass of the system is In the orbital coordinate frame, these vectors with the following expressions The velocity of the point mass m j (j = 1, 2, 3) can be expressed aṡ in which • R j is the derivative of R j with respect to time in the orbital coordinate frame S c .
which is the angular velocity of the orbital coordinate frame, and ν is the true anomaly. Here and below, the dot denotes the derivative with respect to time t. The velocity of particle m j (j = 1, 2, 3) is projected to the orbital coordinate frame and can be expressed as followṡ The kinetic energy of the system is Substituting Eqs. (16), (17), and (18) into the Eq. (19), we obtain the expression of the kinetic energy as The potential energy of the system can be determined as in which R j is the modulus of vector R j , µ is the coefficient of gravitation. Substituting Eqs. (11), (12) and (13) into the Eq. (21) and ignoring the higher order nonlinear terms, the gravitational potential energy can be written as Considering the kinetic energy (20) and potential energy (22), the equations of motion of the three-body tethered satellite system with respect to its center of mass are derived as follows by using Lagrange's equations 3 Relative equilibrium in an elliptical orbit The equations of motions (23) and (24) indicate that there are not any trivial solutions to the equations of motion due to the periodic excitation on the right-hand side of equations. It means the relative equilibrium is impossible for the three-body tethered satellite system in an elliptical orbit. However, it is possible to eliminate the terms on the right-hand side of Eqs. (23) and (24) by regulating the tether lengths L 1 and L 2 , in this case the equations of motion (23) and (24) are homogeneous. As a result, there is a relative equilibrium ϕ 1 = ϕ 2 = 0. The necessary and sufficient condition for the relative equilibrium ϕ 1 = ϕ 2 = 0 is the inhomogeneous terms of equations (23) and (24) are both eliminated. In other words, the tether lengths L 1 (t) and L 2 (t) meet the following differential equations: Let introduce new variables u 1 and u 2 : Then the differential equations (25) and (26) can be written as follows Based on the following equality [37] Eq. (28) be transformed into the following differential equation that the true anomaly is the independent variable: du j dν = e sin ν 1 + e cos ν u j , (j = 1, 2).
The general solution of Eq. (30) as follows in which e is the orbital eccentricity, C j is the integration constant. We assume that the initial value of the true anomaly is ν = 0. Then, considering Eq. (27), the integral constants can be chosen as in which L 1 (0) and L 2 (0) represent the initial value of the tether L 1 and L 2 , respectively. Consequently, if the length of tethers L 1 and L 2 meet the following equations an Earth-pointing relative equilibrium of the three-body tethered satellite system will be obtained for an elliptical orbit with arbitrary eccentricity. According to Eqs. (31), (32) and (33), it is not difficult to verify that there are the following two properties for the tether length control laws: dL j dν = e sin ν 1 + e cos ν L j , (j = 1, 2) .

Hamiltonian of the controlled system
For convenience, new generalized coordinates are introduced as follows: Substituting Eqs. (34) and (36) into Eqs. (20) and (22), the kinetic energy and potential energy of the controlled system can be rewritten as and in which δ is a dimensionless parameter defined as If the generalized coordinates θ 1 , θ 2 and corresponding generalized momenta p θ1 , p θ2 are regarded as canonical variables, the Hamiltonian of the controlled system can be expressed as in which Let us introduce perturbation variables q j , p j (j = 1, 2) using formulas θ j = q j , p θj =νmL 2 2 p j , (j = 1, 2) .
Considering Eq. (29) and the following formula [37] µ r 3 =ν we obtain the Hamiltonian of perturbed motion corresponding to the relative equilibrium ϕ 1 = ϕ 2 = 0, which can be expressed as the following series expansion: in which and the third-order nonlinear term H 3 is eliminated. The coefficients of Eqs. (48) and (49) can be represented as follows

Calculation of the Hamiltonian normal form
Coefficients of Hamiltonian (47) are 2π periodic functions of ν. As a result, the stability of relative equilibrium ϕ 1 = ϕ 2 = 0 can be regarded as the stability of the trivial solution of the nonlinear periodic Hamiltonian system with a period of 2π. According to the stability theorems of periodic Hamiltonian system, the coefficients of the normalized Hamiltonian should be calculated first. An algorithm used to normalize the periodic Hamiltonian with two degrees of freedom is given in [34], using a method for constructing a symplectic mapping. Based on this method, coefficients of the normalized Hamiltonian can be expressed as explicit functions of coefficients of the constructed symplectic mapping. The symplectic mapping is a composite mapping of a linear mapping and a nearly identity mapping. According to the comprehensive description of [34], the linear mapping and the nearly identity mapping, respectively, as follows and q j =q In which, the matrix X (2π) is the value of the fundamental matrix X (ν) of the linearized system at ν = 2π, which elements meet the following differential equations with the initial condition X (0) = E 4 , where E 4 is a fourth-order identity matrix. The generating function S k of the nearly identity mapping (51) can be expressed as in which k = 3, 4 and s i1i2j1j2 = ϕ i1i2j1j2 (2π). Since H 3 in the series expansion of the Hamiltonian (47) is eliminated, the function ϕ i1i2j1j2 (ν) meet ordinary differential equations with initial conditions ϕ i1i2j1j2 (0) = 0. The function g i1i2j1j2 on the right-hand side of Eq. (55) is the coefficient of G 4 , and G 4 is a new Hamiltonian transformed from Eq. (49), using the following univalent canonical transformation in which X (ν) is the fundamental matrix of the linearized system (52). By integrating the ordinary differential equations (53) and (55) numerically in the interval [0, 2π], we can obtain the coefficients of the generating function S 4 .
It is worth pointing out that, in the series expansion of Hamiltonian (46) H 3 = 0, therefore S 3 = 0 in the generating function (54).
It is not difficult to construct a fourth-order symplectic matrix N using the method introduced in [34]. Let replace the canonical variables q j , p j with Q j , P j (j = 1, 2) by the univalent canonical transformation and denoted the transformed generating function as F 4 : According to the conclusions of reference [34], coefficients of the normalized Hamiltonian can be expressed as explicit functions of f i1i2j1j2 which is the coefficient of the generating function F 4 . The normalized Hamiltonian can be represented as follows [34]: (1) If there are no resonances up to fourth-order inclusive, the normal form of the Hamiltonian (47) can be written as (2) If there is a fourth-order resonance, the normal form of Hamiltonian (47) can be written as H = 2 j=1 σ j r j + c 20 r 1 2 + c 11 r 1 r 2 + c 02 r 2 2 + r 1 k 1 2 r 2 k 2 2 (α k1k200 sin γ + β k1k200 cos γ) in which γ = k 1 ϕ 1 + k 2 ϕ 2 − nν; r j and ϕ j (j = 1, 2) are the action-angle variables. The coefficients of the normalized Hamiltonian (58) and (59) are as follows, which are the explicit functions of f i1i2j1j2 in Eq. (57).
On account of H 3 = 0 in the series expansion of the Hamiltonian (47), third-order nonlinear terms are eliminated in the normalized Hamiltonian regardless of whether there is a third-order resonance. Consequently, the stability of third-order resonant case is based on terms of order higher than four in the series expansion of Hamiltonian (47), and an additional approach is required.
The stability analysis of the relative equilibrium begins with the investigation of the linearized system (52). Based on Floquet-Lyapunov theorem [38], in the plane of the dimensionless parameter δ and orbital eccentricity e, the linear stability region and Lyapunov instability region of relative equilibrium are obtained, respectively, in the above three cases. As shown in Figs. 2, 3, and 4, unstable points are distributed along a curve with a maximal or minimal value when the dimensionless parameter δ in the interval (0, 6). Figure 2 and 4 present the stability analysis of relative equilibrium in case 1 and case 3, the unstable curves reach maximal values e = 0.8 at δ = 0.65 and δ = 1.35, respectively, and with the increase of δ the Lyapunov instability curves tend to the line e = 0.6. However, in case 2, the unstable curve reach minimal values e = 0.23 at δ = 1 and with the increase of δ the Lyapunov instability curve tend to the line e = 0.55. Furthermore, in all the above three cases, several Lyapunov instability points occur in the δ − e plane when the orbital eccentricity is close to 1. Obviously, the maximal or minimal value of the instability curve in the δ − e plane is related to the position and mass of the mother satellite and subsatellites, and these values occur in the neighborhood of δ = 1.
There is no asymptotic stability for Hamiltonian systems [35], as a result a further nonlinear stability analysis in the linear stability region is necessary because the stability of the relative equilibrium of the original nonlinear system cannot be derived from the results of stability analysis of the above linearized system. For the stability analysis of nonlinear Hamiltonian system, it should distinguish resonant and nonresonant cases [39].
Based on Arnold theorem [36,39], it can be concluded that if there are no resonances up to fourth-order inclusive and coefficients of the normalized Hamiltonian (58) meet the following inequality then the relative equilibrium is stable. For two degrees of freedom Hamiltonian system, there are at most following 5 types of fourth-order resonance [34]: 4σ 1 = n, 4σ 2 = n, 2σ 1 + 2σ 2 = n, σ 1 + 3σ 2 = n, 3σ 1 + σ 2 = n.
However, there are only two types of fourth-order resonances 4σ 1 = n and 4σ 2 = n for the three-body tethered satellite system in this paper. According to Markeev's theorem [33], in the fourth-order resonant case if coefficients of the normalized Hamiltonian (58) meet the following inequality in which then the relative equilibrium is stable in third-order approximation, otherwise it is Lyapunov unstable. The results of stability analysis of case 1, case 2, and case 3 are presented respectively in Figs. 5, 6, and 7, where the region I represents the stable region of relative equilibrium in the nonresonant case, and II represents the region requiring further analysis. The results of stability analysis indicate that the stable region of the nonresonant case contains all regions of stability in the first-order approximation except the neighborhood of the Lyapunov unstable region of linearized system shown in Figs. 2, 3, 4 and the regions of third-and fourth-order resonance. In the linear stability region, there are large numbers of narrow stability regions close to the line e = 1 and δ = 0; whereas the stability regions are relatively wide in the vicinity of the line e = 0 and the straight line that passes through the point of maximal or minimal value of the unstable curve of the linearized system and is normal to the axis δ. It becomes obvious that among the three cases provided here, the stability region of case 3, shown in Fig. 7, is significantly wider than that of the other two, shown in Fig. 5 and Fig. 6.
In Figs. 8, 9 and 10, we present the results of stability of fourth-order resonance 4σ 1 = n of case 1, case 2, and case 3, where the region I is stable in third-order approximation, the region II represents Lyapunov instability, and III is unstable region of the linearized system. In all three cases investigated here, points of fourthorder resonance 4σ 1 = n are distributed along downward parabolas in the δ − e plane. The vertex of these parabolas lies on a vertical line that passes through the point of maximal or minimal value of the unstable curve of the linearized system and is perpendicular to the axis δ. In the region away from this straight line, the curves of fourth-order resonance 4σ 1 = n are densely distributed. In the above-mentioned three cases, it is obvious that the curves of fourth-order resonance 4σ 1 = n of case 3, shown in Fig. 10, are most sparse, and the unstable points of fourth-order resonance 4σ 1 = n much fewer than that of case 1. As shown in Fig. 8, there are unstable segments in every curve of fourth-order resonance 4σ 1 = n of case 1 in the region e > 0.285 of δ − e plane, while in the region e < 0.285 there is only one unstable point (0.11, 0.57) and the other points of fourth-order resonance 4σ 1 = n are stable in third-order approximation. In the δ − e plane, among the points of fourth-order resonance 4σ 1 = n of case 2 the points (0.51, 0.995), (5.18, 0.88), (5.25, 0.87), (5.28, 0.865) are unstable, whereas the other points of fourth-order resonance 4σ 1 = n all stable in third-order approximation, as shown in Fig. 9. In the δ − e plane, among the points of fourth-order resonance 4σ 1 = n of case 3 the points (0.26,0.925), (3.28, 0.97), (3.48, 0.9258), (3.54, 0.92), (3.6, 0.915) and a segment with end points (4.91, 0.665) and (4.99, 0.635) are unstable, the other points of fourth-order resonance 4σ 1 = n all stable in thirdorder approximation, see Fig. 10.
Results of stability analysis of fourth-order resonance 4σ 2 = n of the case 1, case 2 and case 3 are presented in Figs. 11, 12 and 13, respectively. Different from the distribution of the points of fourth-order resonance 4σ 1 = n, the points of fourth-order resonance 4σ 2 = n only exist in an elliptical orbit with a relatively large eccentricity. In the three cases considered here, only in the region 0.76 < e < 1 of the δ − e plane exists the points of fourth-order resonance 4σ 2 = n and these points are distributed along a curve located in the vicinity of the straight line e = 0.85. Similar to the unstable curve of the linearized system, there is a maximal or minimal value for the curve of fourth-order resonance 4σ 2 = n, and the position of the maximal or minimal value corresponds to that of the unstable curve of the linearized system. Additionally, there are several points of fourth-order resonance near the limit case e = 1. On the whole, most of the points of fourth-order resonance 4σ 2 = n are stable in third-order approximation, and a few points are unstable. In case 1, as shown in Fig Fig. 13. Figures 14, 15 and 16 show cases in which the terms of order higher than four in the series expansion of Hamiltonian (47) need to be considered in the three cases investigated here. These cases include the case of third-order resonance and the case that coefficients of the normalized Hamiltonian (58) meet the following equality c 11 2 − 4c 02 c 20 = 0.
In fact, in the three cases considered here, all the coefficients of the normalized Hamiltonian (58) do not fulfill Eq. (64). While there are two kinds of third-order resonances of 3σ 1 = n and 3σ 2 = n. As shown in Fig. 14, 15 and 16, I represents the region of third-order resonance 3σ 1 = n, II denotes the region of third-order resonance 3σ 2 = n, and III is the instability region of linearized system. Similar to fourth-order resonance 4σ 1 = n, the points of third-order resonance 3σ 1 = n are also distributed along parabolas which vertexes are located in the same straight line with vertexes of the parabola of the fourth-order resonance 4σ 1 = n. While the points of third-order resonance 3σ 2 = n are located in the neighborhood of the instability region of the linearized system. As noted in the previous section, owing to the fact that H 3 = 0 in the series expansion of the Hamiltonian (47), terms of order three in the normalized Hamiltonian are eliminated whether there is a third-order resonance or not. As a consequence, the stability of the third-order resonant case depends on the terms of order higher than four in the series expansion of the Hamiltonian (47).

Conclusions
In this work, the relative equilibrium and stability of a three-body tethered satellite system were investigated, which center of mass moves in an elliptical orbit with arbitrary eccentricity in a central Newtonian field. The system was simplified as the classical two-piece dumbbell model and a control law for regulating the length of tethers was derived, which controlled the system in elliptical orbit to a relative equilibrium. By analyzing the Floquet multipliers of the linearized system, the linear stability region and Lyapunov instability region of relative equilibrium were obtained in the dimensionless parameter plane. In conformity with a method based on the normal form technique and KAM theory, a further nonlinear stability analysis was presented in the linear stability region of the relative equilibrium. Particularly, the stability of relative equilibrium of the three-body tethered satellite system in three cases, including the mother satellite is in the lowest orbit, in the intermediate orbit and in the highest orbit, were investigated in detail.
The results indicate that the relative equilibrium is stable in the linear stability region if there are no resonance up to and including fourth-order. There are two kinds of fourth-order resonances 4σ 1 = n and 4σ 2 = n. In these fourth-order resonance cases, most of the points are stable in third-order approximation, whereas there are only a few Lyapunov instability points in the dimensionless parameter plane. Additionally, there are two kinds of third-order resonances 3σ 1 = n and 3σ 2 = n. However, there are no third-order terms in the series expansion of the Hamiltonian, as a consequence the stability of relative equilibrium in third-order resonant case depends on the terms of order higher than four in the series expansion of Hamiltonian. Comparing the stability charts of the system in the three cases considered here, it is found that in the case of the mother satellite is in the highest orbit and intermediate orbit, the instability points of relative equilibrium in the parameter plane are significantly fewer than that in the case of the mother satellite is in the lowest orbit. In the case of the mother satellite being in the highest or lowest orbit, there are some unstable relative equilibrium points in the region of eccentricity greater than 0.55 in the dimensionless parameter plane. Whereas, in the case of mother satellite is in the intermediate orbit, several unstable relative equilibrium points occur in the region of the eccentricity which is between 0.23 and 0.6 or close to 1 in the dimensionless parameter plane.       6 Results of stability analysis when there are no resonances up to fourth-order inclusive in case 2, in which I is the stable region of the nonresonant case, II is the region requiring further analysis. Fig. 7 Results of stability analysis when there are no resonances up to fourth-order inclusive in case 3, in which I is the stable region of the nonresonant case, II is the region requiring further analysis. Fig. 8 Results of stability analysis of fourth-order resonance 4σ 1 = n in case 1, in which I is the region of stable in thirdorder approximation, II is the unstable region of fourth-order resonance 4σ 1 = n and III is the unstable region of the linearized system. Fig. 9 Results of stability analysis of fourth-order resonance 4σ 1 = n in case 2, in which I is the region of stable in thirdorder approximation, II is the unstable region of fourth-order resonance 4σ 1 = n and III is the unstable region of the linearized system. Fig. 10 Results of stability analysis of fourth-order resonance 4σ 1 = n in case 3, in which I is the region of stable in third-order approximation, II is the unstable region of fourthorder resonance 4σ 1 = n and III is the unstable region of the linearized system. Fig. 11 Results of stability analysis of fourth-order resonance 4σ 2 = n in case 1, in which I is the region of stable in third-order approximation, II is the unstable region of fourthorder resonance 4σ 2 = n and III is the unstable region of the linearized system. Fig. 12 Results of stability analysis of fourth-order resonance 4σ 2 = n in case 2, in which I is the region of stable in third-order approximation, II is the unstable region of fourthorder resonance 4σ 2 = n and III is the unstable region of the linearized system. Fig. 13 Results of stability analysis of fourth-order resonance 4σ 2 = n in case 3, in which I is the region of stable in third-order approximation, II is the unstable region of fourthorder resonance 4σ 2 = n and III is the unstable region of the linearized system. 14 Special points in case 1, where I is the region of third-order resonance 3σ 1 = n, II is the region of third-order resonance 3σ 2 = n and III is the instability region of linearized system. Fig. 15 Special points in case 2, where I is the region of third-order resonance 3σ 1 = n, II is the region of third-order resonance 3σ 2 = n and III is the instability region of linearized system. Fig. 16 Special points in case 3, where I is the region of third-order resonance 3σ 1 = n, II is the region of third-order resonance 3σ 2 = n and III is the instability region of linearized system.