Dynamic behavior analysis of tethered satellite system based on Floquet theory

In this paper, an n-star general dynamic model of tethered satellite system with closed-loop configuration is provided. An analytical method for periodic solution stability of the general dynamic model is proposed based on Floquet theory, which proved that the periodic solution stability of the system depends on the maximum modulus for the eigenvalue of a matrix related to the Jacobian matrix. The periodic solution stability of a 3-star system with equilateral triangle as the initial configuration is analyzed as an example based upon the analytical method, and the results are verified by numerical simulation. The critical spin angular velocity caused by the tether mass and the parameter variation of the 3-star system is analyzed. The results show that the analytical method of periodic solution stability can solve the critical stable spin angular velocity of the tethered satellite system accurately, and the 3-star system can guarantee stable spin in the case of the spin angular velocity is higher than the critical spin angular velocity.


Introduction
Study on tethered satellite system (TSS) has become one central issue of concerns in spacecraft engineering, attracting the attention of researchers in a series of areas [1,2]. Since the first proposal of TSS, it has been extensively applied in actual engineering, such as tethered space robot [3], tethered space harpoon [4], and tethered space net [5]. The first tethered satellite system mission TSS-1 was jointly proposed by the USA and Italy in 1992 [6]. Since then, many research institutions have carried out a large number of space experiments on tethered satellite systems [7] and made abundant achievements [8][9][10][11], but also caused several mission failures [6,12,13].
Based upon this, most researchers concentrated on mechanism design and accordingly the dynamics and control of space capture [14]. An analytical tether length rate control law for the deployment of a flexible TSS under the space perturbations was proposed [15]. Reference [16] proposed a novel adaptive dynamic sliding mode scheme for the deployment of the tethered satellite system. Although the input and command sig-nals are limited based on the strictly bounded terms, the control tension force is discontinuous. Another controller was designed by utilizing small-gain theorem in Ref. [17], while the transversal thrust was considered to suppress the tethers' oscillation. A feedback tension control law for the deployment of space tether with tension constraint and saturation function was proposed by in Ref. [18]. Sakawa [19] used the optimal control to solve the linearized nonlinear system. An input shaping scheme was proposed to generate reasonable planning of the control signal [20], and it is widely used to solve underactuated problem. An earthoriented double-pyramid formation was investigated in Ref. [21] and discovered that the stable configuration was along the nadir direction, and the optimal control of a double-pyramid Earth-facing formation was discussed. A novel adaptive dynamic sliding mode scheme for the deployment of the TSS was proposed in Ref. [15]. Because of the Coriolis force and orbital motion, the research on the spin stability of TSS is essential for the realization of a fast and stable deployment/retrieval [22].
Stability analysis is a core research content in mechanism study of vibration behaviors of dynamic systems. Many stability analysis methods have been used to study nonlinear vibration behaviors of TSS. According to the variation demands of space missions, different satellite formations are often required to complete the missions together. For example, multi-point simultaneous measurements, as required in magnetospheric studies, can be conducted by a set of probes connected by tethers aligned along the local vertical. The potential application of multi-tethers systems is probably greatest in space interferometry both LEO and deep space [23]. Therefore, the study of n-star (or n-body) TSS is of great significance. Correa et al. [24] applied the one-dimensional stable configuration to the 4-body tether system, and pointed out that it can be extended to the tether system of N satellites. Pizarro-Chong et al. [25] studied the tethered formations of hub-and-spoke and closed-hub-and-spoke composed of N bodies, and simulated the system composed of 3 ∼ 6 satellites. The results show that the system composed of more than 4 satellites needs to provide initial spin angular velocity and tether tension to maintain configuration stability. Tragesser [26] proved that the n-star system with closed-loop configuration behaves very similar to a rigid body in the cylindrical equilibrium case when the spin rapid enough, based upon abundant simulation analysis. But there is still a lack of elaboration on the critical angular velocity of the stable spin for TSS with close-loop configuration.
The Floquet theory [27,28] was used to solve a series of dynamic problems [29,30], such as rotor dynamics [31,32], while few researchers used it to study the stability and dynamic behaviors of TSS.
The motivation of this paper is to propose a method for the stable spin conditions of TSS with arbitrary parameter feature based upon a general model with closed-loop configuration. In Sect. 2, a general model of TSS with closed-loop configuration is established and a brief introduction to Floquet theory is also discussed. Section 3 proposes a general method for stability analysis based upon the Floquet theory of the TSS and discusses the numerical results of some examples based on 3-star system. Finally, conclusions and outlooks are drawn in Sect. 4.

General tethered satellite system
A general model of a TSS with n satellites S 1 , S 2 , . . ., S n connected in a closed-loop configuration is considered. The mathematical model of the TSS can be approximated to the bead-spring-damping ring (BSDR) form shown in Fig. 1, no matter the masses of tethers are considered or not, by the following fixed steps. 1. Define m i as the mass of the ith satellite S i ; 2. The ith mass tether connecting satellite S i and S i+1 is divided into finite segment (n i ) units consisting of concentrated mass, massless spring and damping, in which 1 ≤ i < n; 3. The mass m i, j of each unit at the ith segment is defined in which 1 ≤ j ≤ n i ; and spring stiffness k i, j , damping coefficient C i, j , and equilibrium length of each spring L i, j are defined in which 0 ≤ j ≤ n i , based upon the actual rope design parameters.
In order to analyze the dynamic behavior of the system in space, the inertial coordinate system O − XY Z is established by connecting to the centroid of the earth fixedly, and the rotating coordinate system o − x yz is established by connected to the centroid of the tethered satellite system fixedly, respectively, shown in Fig. 2. The positive direction of X -axis is defined as the direction along the centroid of the earth to the ascending  node; the Z -axis is perpendicular to the orbital plane of the TSS system with the positive direction of pointing to the north pole from the centroid of the earth; and the positive direction of Y -axis is determined by the righthand rule, in the inertial coordinate system O − XY Z. The positive directions of y-axis and z-axis are defined as the direction along the centroid of the earth to the centroid of the TSS and the same direction of the Zaxis, respectively, and the positive direction of x-axis is determined by the right-hand rule, in the rotating coordinate system o−x yz. Assuming the revolution angular velocity of the o − x yz around the O − XY Z is , and the positive direction of rotation is counterclockwise.
The motion equation of the ith bead of the BSDR system in the coordinate o − x yz can be expressed as follows in which 1 ≤ i ≤ n +m (including all the lumped mass of n satellites and m = n−1 i=1 n i tether units), and for the ith bread, m i : the mass, r i = (x i , y i , z i ) T : the position vector in the o − x yz, G i : the gravity, T i : the rope tension, D i : the damping force, F e i : the convected inertial force, F k i : the Coriolis inertial force,ṙ i andr i : the first and second order derivative of the position vector r i to the time t, respectively, without considering the perturbation factors such as J 2 perturbation, thermal effect, and atmospheric resistance.
Firstly, the gravitational term in Eq. (1) can be obtained as follows where μ E ≈ 3.986 × 10 5 km 3 /s 2 is the gravity parameter of the earth, r E = (0, −R, 0) T is the centroid position vector of the earth in o − x yz, R = H + R E , in which H is the orbital height of the BSDR system, and R E = 6378 km is the average radius of the earth. Next, the tension term in Eq. (1) , which is determined by the sum of two adjacent tether tensions, can be obtained as follows where r i, j = |r i − r j | is the instantaneous distance between the ith and jth beads, and s i, j satisfies the following relation which control the effective of the tether tension. The damping term in Eq. (1), which indicates the sum of two tether damping along each axis direction based on the linear damping model, can be obtained as the following expression whereṙ i, j is the instantaneous relative speed between the ith and jth beads.
The convected inertial force term in Eq. (1) can be considered to be consisted with the translation convected inertial force and the centrifugal convected inertial force, which can be expressed as follows where is the value of spin angular velocity of the coordinate o − x yz which equals to its revolution angular velocity numerically, and ω = (0, 0, ) T andω represent the rotational angular velocity and angular acceleration of the o − x yz, respectively.ω = 0 must be hold by supposing the orbit of revolution is circular, then the relation can be obtained as follows where i, j , and k represent the x, y, and z axes unit vectors in the o − x yz. Furthermore, the convected inertial force can be repressed as Finally, consider the Coriolis inertia force term in Eq. (1) can be obtained as whereω = ω is the spin angular velocity of the o−x yz, then Eq. (10) can be rewritten as

Brief introduction to Floquet theory
Floquet theory is a stability theory for a periodic variable coefficient linear ordinary differential equation, which was put forward by G. Floquet in 1883 [33]. This section will discuss the general method for N dimensions equation. Consider a system of first-order differential equationṡ where A is a matrix with period T . Supposing the system (12) has N linearly independent particular solutions and the initial conditions are satisfied It is easy to prove that this set of particular solution is the fundamental solution system of system (12). Based another particular solution, which can be expressed as in which (1 ≤ i, j ≤ N ). Substituting t = 0 and Eq. (14) into Eq. (15), the equation can be solved as follows In order to determine the stability of the periodic solution of the system (12), it is necessary to find a particular solution which satisfies the following conditions in which p * (t) and p * (t + T ) can be expressed as the fundamental solution system (12) and respectively. It can be obtained by utilizing Eqs. (15), (18), and (19) into Eq. (17) as follows which is equivalent to by assuming and It means that the following relation must be held Ulteriorly, Eq. (24) can be rewritten as follows by assuming and where ρ is the eigenvalue of the matrix A T , which means that the zero solution is asymptotically stable, if and only if all |ρ| < 1, otherwise the zero solution is unstable.

General method for stability analysis of TSS
The motion Eq. (1) is rewritten into a first-order form, to analyze the periodic stability of TSṠ where and Assuming P 0 : (x) is a fixed point on the Poincare section, wherē . . , p 6(n+m) T consisted of 2 (n + m) small amounts, and P 1 : (x + p) is a point in the small neighborhood of P 0 : (x) that P 1 can be regarded as a perturbation of P 0 . Therefore, P 1 can be expanded into Taylor series form by retaining the linear terms based upon P 0 , which can be obtained as followṡ where The system (36) with the initial conditions can be solved by dividing into the following two cases, based on the ordinary differential equation theory.
1. When the matrix J has N = 6 (n + m) linearly independent eigenvectors v 1 , v 2 , . . . , v N , and their corresponding eigenvalues are λ 1 , λ 2 , . . . , λ N . Then a basis solution matrix of homogeneous linear differential equations with constant coefficients Eq. (28) can always be expressed as the general solution can be expressed as where C = (C 1 , C 2 , . . . , C N ) T is a matrix consisted of undetermined constants. Substituting the initial condition (36) into Eq. (38), the equation can be obtained as follows is the eigen-matrix of J | x=x , and the constant matrix satisfies the relation as follows Thus, the general solution can be rewritten as 2. When the matrix J has less than N linearly independent eigenvectors, it is easy to prove the following conclusion. If λ i is the N i -fold eigenvalue of J, then there must be N i linearly independent solutions in the form of in which v i,0 is a non-zero solution of and a recurrence relation can be obtained as follows where 1 ≤ j ≤ N i − 1. Define the matrices P and V are separately consisted of N linearly independent vectors p i and v i , then the general solution of the system in this situation can be solved by substituting P and V into (42).
The eigenvalues of the matrix p (t) at the periodic sampling time ρ 1 , ρ 2 , . . . , ρ N are the Floquet multiples, according to the Floquet theory in Subsect. 2.2. It is easy to prove that the Floquet multiplier of the system is the eigenvalue of the matrix by assuming t 0 = 0, and the periodic solution must be asymptotically stable if and only if the following condition is satisfied 3.2 Verification based on a typical example Taking a 3-star system connected by massless tethers as an example, the stability analysis is carried out based on the above general model and analysis method, by letting The motion equation of the system can be obtained as follows where m i = M, E i, j = E, A i, j = A, C i, j = C, L i, j = L, and i = 1, 2, 3, by assuming the initial position and velocity are in the xoy plane. Considering that the system has an initial steady-state configuration of regular N -polygon, which may be assumed that the 1st satellite is located on the y-axis at the initial moment, the distance from the origin of the o − x yz is R 1 . Then the initial position of any satellite can be obtained as follows in which N is the number of satellites, and δα = 2π N expresses the angle between two adjacent satellites. It is easy to prove that R 1 satisfies the following relation in which ω expresses the spin angular velocity of the TSS system. In order to verify the initial steady-state configuration of tethered satellite system, a simulation model was established. The trajectories of the 1st, 2nd, and 3rd satellites are shown in Figs. 3 and 4, which are represented by red, green, and blue dotted curves, respectively; and the black circles and the solid lines in black represent the satellites and the tethers connecting two adjacent satellites, under the spin angular velocity condition |ω| = 2.0 and |ω| = 2.5 , separately. Based on the motion Eq. (49) and the initial conditions (50), the TSS can maintain the initial configuration depend on its spin under both conditions, which proves the correctness of the steady-state configuration and the solution of the initial conditions.
The curve of Floquet multiplier depended on clockwise spin angular velocity of the TSS can be obtained based upon the method mentioned in Subsect. 3.1 with the parameters shown in Table 1. It can be found from Fig. 5 that the configuration is stable when the spin angular velocity is higher than the critical  Obviously, the configuration of the disturbed system cannot be held when the spin angular velocity is lower than the critical angular velocity (| ω |= 2.0 < 2.17), while the disturbed system will converge to the original configuration in about 2 periods, when the spin angular velocity satisfies | ω |= 2.5 > 2.17.
The relative deviation from the initial steady-state configuration is defined as follows to analyze the convergence of each satellite to the initial configuration depends on time. The relative deviation curves of two operating modes with the two spin angular velocities of |ω|= 2.0 and |ω|= 2.5 are shown in Fig. 8. It is shown that the relative deviation of each satellite in the disturbed system cannot be attenuated to 0 when the spin angular velocity is less than critical angular velocity, shown in Fig. 8a; but the relative deviation of each satellite will gradually decay to 0 to restore the initial steady state configuration in a limited time, when the spin angular velocity satisfies | ω |= 2.5 > 2.17. Therefore, the solution of periodical stability obtained by the general method based upon the Floquet theory mentioned in Subsect. 3.1 has been proved that the configuration solved by (50) and (51) is stable when the clockwise spin angular velocity rapider than about 2.17 times the revolution angular velocity.

Sensitivity analysis of system parameters
The effects of the critical angular velocity caused by the variation of the system parameters are considered based upon the above 3-star system. Firstly, two groups of 3-star system with different satellite masses are analyzed, in which the satellite masses in the first group are 50 kg (less than the original system) and which in the other group are 1000 kg (greater than original system), respectively. It can be found from Fig. 9 that the changing of satellite mass will significantly affect the critical angular velocity of the stable spin of the system, based upon the Floquet theory. The Floquet multiplier of the system with the satellite masses of 50 kg and 1000 kg are shown in Fig. 9a, b, c, and d, respectively, and the critical angular velocities are separately about 2.79 and 2.07 times of the revolution angular velocity. This leads a conclusion that the critical angular velocity will decrease with the increase in satellite mass, which means that the increase in the system inertia caused by increasing satellite mass will reduce the stable spin condition. The critical angular velocities solved by the Floquet theory are verified by disturbed simulations. The configuration of the disturbed system with the satellite mass of 50 kg, shown in Fig. 10, can hold the original configuration in about 2 periods, if the spin angular velocity satisfies | ω |= 2.9 > 2.79. As shown in Fig. 11, the configuration of relative deviation of each satellite in the disturbed system will gradually decay to 0 to restore the initial steady state configuration in a limited time when and only when the spin angular velocity is larger than the critical angular velocity, shown in Figs. 10c, f and 11c, f. Secondly, the applicability of the general model and its stability discrimination method for complex system with multi-body is demonstrate. Considering the 3-star system consisted of multi-beads (more than 3) with different masses, introduced by the mass effect of tethers. Assuming the tether connecting two satellites with the mass satisfies M = 1000 kg can be divided into two segments with the center mass is m = 10 kg. Considering that each tether is divided into two segments, and the steady-state configuration of 3-star system with 6 beads is shown in Fig.12. The parameters R 1 and R 2 can be determined by the following equations The critical angular velocity of the 6-bead system can be determined by analyzing the Floquet multiplier in Fig. 13, which implies a conclusion similar to that shown in Fig. 9d, that is |ω|≈ 2.08 . The simulation results of the stability of the disturbed 6-bead system under different spin angular velocities are shown in Fig. 14. It can be found that the simulation results verified the correctness of the stability condition implied by the results of the Floquet multiplier shown in Fig. 13. The trajectory of each bead of the disturbed system will gradually approach a circular orbit in about 10 periods in order to stabilize the system configuration, when the spin angular velocity is larger than the critical angular velocity (|ω|≈ 2.08 ). On the contrary, the system cannot maintain configuration steadily when the spin angular velocity is not large enough. However, the 6-bead system cannot approach the initial stable configuration in a limited time even when the spin angular velocity is greater than the critical angular velocity. It is means that a rapid enough spin can keep the configuration stable but cannot completely suppress the vibration, for the system with more than 3 beads. In other words, it is obvious from Fig. 15a and c that the relative deviation of the center distance of each beads cannot approach to 0, whether or not the spin angular velocity is greater than the critical angular velocity; but the relative deviation of the tether length can approach to 0, only when the spin angular velocity is greater than the critical angu- Fig. 11 The simulation results of the system with the satellite mass of 1000 kg. a for the orbit of each satellite when the clockwise spin angular velocity of the system satisfies |ω|= 1.5 , b for the top view of (a), (c) for the relative deviation of the system sat-isfies |ω|= 1.5 , d for the orbit of each satellite when the clockwise spin angular velocity of the system satisfies |ω|= 2.2 , e for the top view of (d), and f for the relative deviation of the system satisfies |ω|= 2.2

Fig. 12
The steady-state configuration of 3-star system with 6 beads lar velocity. Considering there are multiple solutions of the geometry of the system with the constant tether length but different included angles, when the number of beads in the system exceeds 3. The kinetic energy causing the change of the included angle of adjacent tethers will not be dissipated, since there is no constraint on the rotational relationship between each bead and tether. Therefore, it can be concluded that although the spin stability conditions of the closed-loop system with more than 3-bead can be accurately analyzed by the general model and the method proposed in Fig. 13 The Floquet multiplier of the 6-beads system this paper, the periodic change of configuration caused by the variation of tether angle cannot be suppressed only by accelerating the spin angular velocity.

Conclusion and outlooks
In this paper, a general dynamic model of n-star tethered satellite system with closed-loop configuration has been proposed. An analytical method of periodic solution stability of the general dynamic model is provided based on Floquet theory, which proves that the periodic solution stability of the system depends on the maximum modulus of the eigenvalue for a matrix related to the Jacobian matrix of the system. On the basis of the general periodic solution stability analysis method, taking a 3-star system with the initial configuration of equilateral triangle, the critical stable clockwise spin angular velocity of the system is analyzed, and the numerical simulation is carried out to verify the results. The results show that the analytical method of periodic solution stability can solve the critical stable spin angular velocity of the TSS accurately. The system can guarantee stable spin when the value of its spin angular velocity is higher than the critical stable spin angular velocity; otherwise, the disturbed system will not be able to re-converge to the initial configuration in finite time. The periodic oscillation of the stable configuration consisted of more than 3 beads, caused by variation of the included angle of the tethers, cannot be suppressed only by accelerating the spin angular velocity. Therefore, the vibration suppression scheme except changing the spin angular velocity is of great significance for the multi-bead (more than 3) system, which is being actively studied by the authors.

(a) (b)
(c) (d) Fig. 15 Relative deviation of the system. a for the distance from each bead to the center of the system with the spin angular velocity of |ω|= 1.5 , b for the length of each tether with the spin angular velocity of |ω|= 1.5 , c for the distance from each bead to the center of the system with the spin angular velocity of |ω|= 2.2 , and (b) for the length of each tether with the spin angular velocity of |ω|= 2.2