Dynamics and Control of Tethered Satellite System in Elliptical Orbits under Resonances

School of Mechanical Engineering, Institute of Launch Dynamics, Nanjing University of Science and Technology, Nanjing 210094, China School of Mechanical Engineering, Nanjing University of Science and Technology, Nanjing 210094, China MOE Key Laboratory of Impact and Safety Engineering, Ningbo University, Ningbo 315211, China School of Science, Nanjing University of Science and Technology, Nanjing 210094, China


Introduction
Tethered satellite system is a promising new type of spacecraft [1][2][3]. It has great potential in applications such as space debris capture [4], space elevator [5], and orbital transfer [6]. For different applications, there will be different combinations of system parameters, which will increase the possibility of resonance. Especially for the system in elliptical orbit, the resonance of the system becomes more complex due to the periodic excitation. Moreover, the dynamic characteristics of the resonance system are very different from those of the nonresonance system, so it is necessary to study all kinds of resonance motion of the tethered satellite system.
The researches on the dynamic behavior of the tethered satellite system under nonresonance mainly focus on the steady-state solutions and stabilities. Takeichi et al. [7] studied the periodic solution of the librational motion of a tethered system in elliptic orbit via the Lindstedt perturbation method. Their analysis showed that the periodic solution is the minimum energy solution, and from the mechanical point of view, the periodic solution in an elliptic orbit has the same significance as the equilibrium state in a circular orbit. Sidorenko and Celletti [8] obtained the family of planar periodic motions for the "spring-mass" model of tethered satellite systems and studied the bifurcations and the stability of these periodic motions concerning in-plane and out-ofplane perturbations. Nakanishi et al. [9] investigated the inplane periodic solutions of a dumbbell satellite system in elliptic orbits by bifurcations with respect to the orbital eccentricity. The periodic solutions are projected on the van der Pol plane to observe the directional change points of the trajectories which can become a useful tool to develop new control schemes. Peláez and Lara [10] employed the Poincaré method of continuation of periodic orbits to obtain other periodic solutions in electrodynamic tethers on inclined orbits that cannot be obtained by using asymptotic techniques. Burov et al. [11] studied the families of periodic motions analytically in a cabin-elevator system, and the condition for the existence of periodic solutions was determined.
Usually, the periodic motions of the system are the desired target motions, but most of these periodic motions are unstable; many control methods are proposed to stabilize the unstable periodic motions. Fujii et al. [12] used a timedelay feedback control (DFC) to control chaotic librational motion of a tethered satellite in an elliptic orbit; the result shows that the DFC scheme is particularly effective in controlling the chaotic motion of the system. The same control method is employed by Peláez and Lorenzini [13] to control an electrodynamic tethered system in inclined orbit. But this control law does not stabilize the unstable periodic orbit for reasonable values of the control parameters. Thus, for the librational control of the electrodynamic tethered system, Iñarrea and Peláez [14] adopted the extended time-delay autosynchronization (ETDAS), while Williams [15] used the predictive control with time-delayed feedback scheme; all achieved good control performances. Additional, Kojima et al. [16] applied a model-following, decoupling-control method combining with the delayed feedback control method in a new approach to control the librational motion of the tethered satellite system in elliptic orbits. The numerical simulations showed that the control method has good performance such as a short settling time for convergence and small control forces.
However, there are few studies on the nonlinear resonances of a tethered satellite system, and the limited researches mainly focus on the system in circular orbit [17,18], and the more complex elliptical orbit system which may have combination resonance needs to be studied. This paper studies the nonlinear resonance of TSS in elliptical orbits. The mathematical model of the in-plane tethered satellite system in elliptical orbits in which the main satellite is treated as a rigid body is derived in Section 2. Then, perturbation analysis of the system is executed to get the parametric relations of every kind of resonance in Section 3. In Section 4, the bifurcation analyses of the TSS under resonance are studied numerically. Section 5 implements the ETDAS technique to stabilize the unstable motion under resonance to a periodic motion. Periodic solutions under different resonances are compared in Section 6. The conclusions are given in Section 7.

Mathematical Model
Consider an in-plane tethered satellite system as shown in Figure 1, where the main satellite is treated as a rigid body of mass M and the subsatellite is envisioned to be a point of mass m. The subsatellite is connected with an inelastic massless tether l at a joint point of distance ρ to the mass center C of the main satellite. The mass of the main satellite is assumed to be much greater than the mass of the subsatellite, such that the center of mass of the system can be assumed to coincide with the main satellite that moves in an unperturbed Kepler elliptical orbit of semimajor axis a and eccentricity ratio e.
Two right-oriented reference frames are used to describe the system motion. The first is the earth-centered inertial frame denoted by O-XYZ, the origin of which is located at the center of the Earth, where i, j, and k are unit vectors in the directions of OX, OY, and OZ, respectively. The second is that the body frame established with the directions of the axes coincided with the principal axes of the main satellite. I x , I y , and I z are the principal moments of inertia of the main satellite expressed in the body frame. In this paper, only a specific relation I y = I x is considered.
According to Figure 1, the position vectors of the main satellite and the subsatellite in the inertial frame are, respectively, where ν is the true anomaly and R is a distance between the main satellite and the center of the earth. Because the main satellite moves in an unperturbed Kepler elliptical orbit, it is easy to get R = að1 − e 2 Þ/ð1 + e cos νÞ. The potential energy of the system is in the following form [19]: where μ e is the gravitational constant. The kinetic energy of the system is The equation of motion of the system can be derived by Lagrange's equations, Considering that the tether length is much smaller than the orbit radius, the joint point of distance ρ is smaller than the tether length, and only a small eccentricity ratio e is studied here, the equations of motion of the system are obtained as follows: where β 1 = ρ/l, β 2 = Mρl/I z , and the overdot denotes the derivative with respect to the true anomaly ν.

Perturbation Analysis
The resonance motions of the system are analyzed by a multiscale method. It is assumed that the eccentricity is small, and a small scaling factor ε is introduced to rescale the eccentricity as follows: The first-order approximate analytic solutions of Equation (6) can be expressed as where T i = ε i τ, i = 0, 1, represents different time scales. By substituting Equation (8) into Equation (6) and equating the same power of ε, one can get the following: (i) Order (ε 1 ) (ii) Order (ε 2 ) where D r = ∂/∂T r ðr = 0, 1Þ represents the partial differential operators.
The first-order approximation solutions of Equation (9) take the complex form where j is an imaginary unit and the coefficients A 1 ðT 1 Þ and A 2 ðT 1 Þ are functions of time T 1 . And cc denotes the complex conjugate of the preceding term. The coefficients B 1 , B 2 , and Γ 1 are time-independent, their expressions are as follows: And jω 1 and jω 2 are the two mutually different pure imaginary roots of the following characteristic equation of λ: So it is easy to get ω 2 1 = 3 + 3β 1 and ω 2 2 = 3β 2 . Substituting Equation (11) into Equation (10) yields where the coefficients R i and P i ði = 1, 2, ⋯7Þ are shown in the appendix.

International Journal of Aerospace Engineering
It is easy to see that the exponential terms in the righthand side of Equations (14) and (15) So, it can be seen that internal resonances and parametrically excited resonances may occur in the system. All possible resonance types of the system and their corresponding parameter relationships are given in Table 1. In the table, cases 1-3 are internal resonances and cases 4-8 are parametrically excited resonances.
As can be seen from Table 1, when a certain parametric relationship is satisfied, the system may have the internal resonance or parametrically excited resonance. β 1 is always positive, so one can get ω 1 > ffiffi ffi 3 p . Therefore, case 6 does not exist, that is, combination resonance 3 (CR3) will never occur in the system.
From Table 1, in cases 1-5, the relationship between parameters β 1 and β 2 is positively correlated. Figure 2 shows the relation curves between parameters β 1 and β 2 under every resonance type. Considering that most applications of TSS have ρ < l, the range of β 1 in the figure is ½0, 1. It can be seen from the figure that, except for SR2, that is β 1 = 1/3; the value range of parameter β 2 is always in [0,8]. That means, in most applications of space tethered systems (0 < β 1 < 1), including elliptic orbits, the resonance region can be avoided when β 2 > 8. It can provide a reference for the parameter design of TSS.

Bifurcation Analysis
4.1. Bifurcation Diagram Changing with Orbital Eccentricity e . In Section 3, the parametric relations of every kind of resonance are obtained by the analytic method. In this section, the bifurcation behaviors of the system under nonresonance (NR) are compared with those under resonances by numerical simulations. Consider that β 1 is usually a small quantity; β 1 = 0:001 is chosen in the simulations. So ω 1 ≈ ffiffi ffi 3 p can be obtained. Table 2 shows the values of parameter β 2 corresponding to every resonance type. It can be seen that when β 1 = 0:001, six kinds of resonances can occur in the system by choosing the value of parameter β 2 , where SR2 and CR3 will not occur. And notice that these six kinds of resonance will not occur at the same time.

Case 1. NR
Firstly, the bifurcation behavior of the system under NR with orbital eccentricity e is investigated. The parameters β 1 = 0:001 and β 2 = 10 are selected. No resonance occurs in the system. Figure 3 shows the bifurcation diagrams of the system changing with orbital eccentricity e. As can be seen from the figure, the bifurcation diagrams of θ and α are similar. In e ∈ ½0, 0:284Þ, they experience bounded motions of periods, period-doubling, and quasiperiods. When e ≥ 0:284 , the motions of θ and α are all chaotic. As we all know, for a dumbbell model, the pitch motion θ enters a tumbling state around the eccentricity of 0.313. So when the mother satellite is treated as a rigid body, this system enters a chaotic motion easier than a dumbbell model. But the critical orbital eccentricity does not decrease much. The coupling of the motions of θ and α makes it easier for the system to enter chaotic motions. But the coupling is not very strong in the case of nonresonance.

Case 2. CR1
Next, the bifurcation diagrams of θ and α are observed under combination resonance. A set of parameters is selected as β 1 = 0:001 and β 2 = 0:18; one can get ω 1 ≈ ffiffi ffi The bifurcation diagrams under CR1 are shown in Figure 4. As can be seen from the figure, the bifurcation diagrams of θ and α show significant differences. The bifurcation diagram of pitch angle θ is similar to that in the case of nonresonance, and the critical value of e is 0.279. However, the bifurcation diagram of α is quite different from that of the nonresonant case. The motion of α enters into chaotic motion at a very small value of e, and the critical value of e is just 0.08. And it is worth noting that the amplitudes of the periodic, doubling-periodic, and quasiperiodic motion of α also increase significantly compared with those of the nonresonance. This type of parametrically excited resonance has a great impact on the motion of α, while it has a small impact on the motion of θ.
International Journal of Aerospace Engineering Case 3. PR In this case, the bifurcation diagrams under primary resonance are obtained. Here, β 1 = 0:001 and β 2 = 1 are selected; one can get ω 1 ≈ ffiffi ffi 3 p and ω 2 = ffiffi ffi 3 p , so ω 1 ≈ ω 2 . Figure 5 shows the bifurcation diagrams under the primary resonance. As can be seen from the figure, the bifurcation diagrams of θ and α are similar to those in the case of combination resonance. The critical values of e are 0.282 and 0.133 in the bifurcation diagrams of θ and α, respectively.
After simulating the bifurcation diagrams of the other four types of resonances at β 1 = 0:001, the results similar to the above two kinds of resonances can be obtained. That is, the resonances have a great impact on the motion of α, while they have a small impact on the motion of θ. The main effects of resonance on the motion of α are that the amplitudes of its periodic, doubling-periodic, and quasiperiodic motions increase compared with those of nonresonance. The critical value of orbital eccentricity that causes it to move into a chaotic motion is much smaller than that of nonresonance. The main reason for this situation can be obtained by observing the system's first-order approximate differential Equation (9). In the equation, the motion of θ decouples with α. So internal resonances and parametrically excited resonances related toαhave very little effect on it. The first-order approximation differential equation of α in Equation (9) is coupled with θ, so the motion of α is greatly affected by these kinds of resonances.

Case 4. SR2
By observing Equation (9), it can be seen that the motion of θ is affected by the parametric excitation term, so it can be suspected that the SR2 may have great effects on the motion of θ. A set of parameters is selected as β 1 = 0:33 and β 2 = 10 in bifurcation simulations; one can get ω 1 ≈ 2 and ω 2 = ffiffiffiffiffi 30 p , so the system is under SR2. The bifurcation diagrams under superharmonic resonance are shown in Figure 6. As can be seen from the figure, the bifurcation diagrams of θ and α are similar. The bifurcation diagram of θ is significantly different from those of nonresonance and other types of resonances. The main manifestation is that the amplitudes of the periodic, doubling-periodic, and quasiperiodic motions of θ are significantly larger than those in other cases when the eccentricity is small. Moreover, the critical value of e that makes the motion of θ be chaotic is also reduced to 0.204. As the motion of α is coupled with that of θ, their bifurcation diagrams are similar. The critical value of e is also 0.204 in the bifurcation diagram of α.

Bifurcation Diagram
Changing with Parameter β 2 . The bifurcation behaviors of the system changing with parameter β 2 are investigated next. Given that β 1 = 0:001 and e = 0:15, Figure 7 shows the bifurcation diagrams changing with parameter β 2 . There are six vertical solid lines in the figure; they are the values of parameter β 2 under various kinds of resonance when β 1 = 0:001 as shown in Table 2. It can be seen in the figure that the change of parameter β 2 has little effects on the motion of θ, that is, the six types of resonance listed in Table 2 have little effect on the motion ofθ. The motion of α varies significantly with the change of parameter β 2 , and the regions with a large amplitude of α correspond to the positions of resonance regions given in Table 2. So the resonance regions obtained by analytical analysis are in good consistency with the results obtained by numerical simulations. This is also consistent with the above results of bifurcation behaviors with the change of parameter e. In Figure 7, it can also be seen that there are two areas with an obvious sudden increase of amplitude in the bifurcation diagram of α, which are not corresponding to the above resonance regions. The two areas are located in the center of the region β 2 = 3 and β 2 ≈ 4:64, respectively. They correspond to superharmonic resonance which is ω 2 = 3 and combination resonance which is ω 2 − ω 1 = 2. The two resonances are ignored in the theoretical analysis because the dynamic equations are truncated by order 2 in the previous perturbation analysis.
The above analysis shows that SR2 has a greater effect on the motion of θ. Here β 1 = 1/3 is selected, so ω 1 = 2, and SR2 occurs in the system. Table 3 shows the values of parameter  The value of ω 2

Frequency relationship
Types of resonance 5 International Journal of Aerospace Engineering β 2 corresponding to every resonance type at β 1 = 1/3. It can be seen that when β 2 takes the values in Table 3, multiple resonances occur simultaneously in the system. Figure 8 shows the bifurcation diagrams changing with parameter β 2 at β 1 = 1/3. It can be seen that the bifurcation diagrams of θ and α are similar, and both of them have a sharp increase in amplitude in the resonance regions. It also confirms the above results of bifurcation behaviors with the change of parameter e.
From the above, since the first-order approximate equation of θ is not coupled with α, only SR2 has a great effect on its motion. However, the amplitude of the motion of α is significantly increased under various kinds of resonances and can even be chaotic. In short, resonances have a huge effect on the motion of the system.

Periodic Control
According to the analysis in the previous section, resonances will make the system more prone to chaotic motion, which will bring harm to the operation of the system. In general,  International Journal of Aerospace Engineering we can make the system avoid the resonance region by selecting appropriate system parameters, but sometimes, when the cost of changing system parameters is higher, control schemes are needed to make the system operate normally. Therefore, this section discusses a control strategy to make the chaotic motion of the system to the periodic orbit in the resonance region. The extended time-delay autosynchronization (ETDAS) has been applied successfully to stabilize the chaotic motion of tethered satellite systems. This control technique has some great advantages that make it widely used. Firstly, it only requires the knowledge of the period of the desired periodic orbit instead of a reference signal of the desired regular motion. Secondly, when the system operates in the neighbor-hood of the desired periodic solution, the control inputs will take small values. And this method does not require fast switching or sampling that makes it easy to apply in practice. Based on these excellent advantages, the ETDAS will be employed to stabilize the unstable periodic motions in this paper. Assuming that the system is acted upon by additional forces regardless of the actual actuators, the controlled governing equation becomes where the feedback control inputs F 1 ðνÞ and F 2 ðνÞ are given by   Table 3: The values of parameter β 2 corresponding to every resonance type at β 1 = 1/3.

The value of β 2
The value of ω 2

Frequency relationship
Types of resonance 7 International Journal of Aerospace Engineering where k 1 and k 2 are the feedback gains and R 1 and R 2 are the memory parameters. The delayed time τ is periodic of the desired periodic orbit.
The basic block chart of the ETDAS control scheme is shown in Figure 9.
Pang and Jin [20] have verified the validity of ETDAS for the TSS studied in this paper when it is in a nonresonant region by numerical and experimental methods. Therefore, this section mainly investigates whether the ETDAS method is still valid in resonance regions. The stability analysis method proposed by Bleich and Socolar [21] is adopted, the stability domains of the controlled tethered satellite system have been calculated. For simplicity, the study of the stability is limited to the case of R 1 = R 2 = R and k 1 = k 2 = k. As previous studies have shown, the system becomes harder to control as the orbital eccentricity increases. Therefore, there is a threshold value e * of orbital eccentricity. When the threshold value is exceeded, the method cannot control the system to a desired periodic orbit. The stability domains will be shown as functions of the control parameters R, k, and e * . Without loss of generality, three cases analyzed in the previous section are selected for comparison, namely, NR, CR1, and SR2. Figure 10 gives the stability domains provided by the ETDAS method in the three-dimensional space (R, k, and e * ) for the three cases. In the figure, the upper bound surfaces of the stability domains are given for various cases, and the stability domains are the regions between each surface and the (R, k) plane of e * = 0. The top surface is NR, the middle one is CR1, and the bottom one is SR2. It can be seen from the figure that the stability domains of resonance cases are smaller than that of nonresonance one. Among them, the largest threshold value of e is e * = 0:47 for NR and e * = 0:41 for CR1, and the smallest one is e * = 0:09 for SR2. Besides, it can be seen from the figure that the control stability domains are not particularly sensitive to the change of parameter R. Therefore, in order to show the changes of the stable domains more clearly, Figure 11 shows the stable domains in two-dimensional space (k, e * ) when R = 0:5. As can be seen from the figure, the stable domain becomes smaller in resonance cases than in nonresonance cases, among which the SR2 is reduced very much. It can also be seen that when the control parameter k is small, it is of great help to increase the threshold value of eccentricity. When k is large, its increase can hardly increase e * . So the control purposes can be achieved without setting the value of k too large.

Periodic Solutions under Resonances
The multiscale method described in Section 3 of this paper can be used to obtain the periodic solution of the system in the case of small orbital eccentricity. But for every case of resonances, the method has to be recalculated to find every analytic solution. To simplify the analysis, a straightforward expansion method is used to obtain a 2π basic periodic solution in this section. This method has been successfully used by Peláez et al. [22] to obtain the 2π periodic solution of an electrodynamic tethered system. The basic periodic solutions of the system studied in this paper have been given by Pang et al. [23]. The correctness of the analytical solutions in the nonresonant case has been verified. To prove that the  Figure 8: Bifurcation diagrams changing with parameter β 2 at β 1 = 1/3.

Dynamical system
Delays j Input Output Figure 9: Basic block chart of the ETDAS control method. 8 International Journal of Aerospace Engineering solutions are still reliable in resonance, the numerical and analytical solutions of CR1 and SR2 are compared. Figure 12 shows the time histories of θ and α in two orbits corresponding to different orbital eccentricities under CR1. Three cases are given in the figure corresponding to the following values of the orbital eccentricity: e = 0:01, 0:03, 0:05. In the figure, the curve with a larger amplitude corresponds to higher orbital eccentricity, and the straight line corresponds to the numerical solutions and the circle to the analytical solutions. It can be seen from Figure 12 that, within two orbital periods, the analytical solutions of periodic motion are in good agreement with the numerical solutions, and the smaller the eccentricity value is, the better the coincidence is. This shows that the approximate periodic solutions obtained by the perturbation method are correct and reliable in this kind of resonance.
Then, the time histories of θ and α in two orbits corresponding to different orbital eccentricities under SR2 are shown in Figure 13. The meanings of the curves in the figure are consistent with those in Figure 12. It can be seen from Figure 13 that, within two orbital periods, the analytical solutions of periodic motion are in good agreement with the numerical solutions, indicating that the approximate periodic solutions obtained by the perturbation method are correct and reliable under such resonance conditions. It is worth noting that the amplitudes of the periodic solutions corresponding to different orbital eccentricity are larger than those under CR1.
To observe the differences in periodic solutions under different resonance conditions, Figure 14 shows the phase diagram of basic periodic solutions in phase space in the same orbital eccentricity e = 0:05 under NR, CR1, and SR2. It can be seen from the figure that SR2 has the largest amplitude of periodic solution, CR1 has the second-largest amplitude, and NR has the smallest amplitude. The amplitude of the periodic solution in SR2 is about ð0:4, 0:46Þ, CR1 is about ð0:05, 0:28Þ, and NR is about ð0:05, 0:06Þ. The amplitude is much larger in resonance than in nonresonance. To further illustrate, Figure 15 shows the phase diagram of the basic periodic solution in phase space under three cases of NR, CR1, and SR2 at the same orbital eccentricity e = 0:1. As can be seen from the figure, the amplitude of SR2 is about ð 0:71, 1:45Þ, CR1 is about ð0:1, 0:56Þ, and NR is about ð0:11 , 0:12Þ. Obviously, SR2 has a much larger amplitude than NR, especially the motion of α.
The results from the above phase diagram simulations show the major cause for the decrease of the control ability of ETDAS under resonance, because the purpose of ETDAS is to control the motion of the system to a periodic orbit. As can be seen from the above phase diagram simulations, the amplitudes of periodic motion of the system in resonance are larger than those in nonresonance. Even in some resonance cases, the amplitudes increase very much, which directly makes the system more difficult to control in resonance than in nonresonance.

Conclusion
In this study, resonance motions of tethered satellite systems in elliptical orbits are studied. The perturbation analysis of the system is carried out by using the multiscale method. All possible resonance types and corresponding parameter relations of the system are given, including internal resonances and combination resonances. The resonance parameter domain which is in 0 < β 2 ≤ 8 is given to provide a reference for parameter design of a space tethered system. The bifurcation behaviors of the system with orbital eccentricity e and parameter β 2 are studied by the numerical method. The results show that the threshold values of e in resonance cases are much smaller than those in nonresonance cases, that is to say, it is easier to enter into chaotic motion under resonances. And the resonance regions    The ETDAS method is applied to stabilize the chaotic motion of the system. Stability analysis shows that stable domains become smaller in resonance cases than in nonresonance cases. The 2π periodic solutions are obtained analytically. By comparing the phase diagrams, it can be observed that the amplitudes of the periodic solutions in resonances are much larger than those in nonresonances. This is the main reason for the decrease of the control ability of ETDAS under resonances. In future work, some new control strategies should be designed so that the system can still be controlled to periodic motion from a chaotic one in high orbital eccentricity under resonances.

Data Availability
Data will be available on reasonable request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.