# Dynamics and Control of a Tethered Satellite System Based on the SDRE Method

Authors: Yong-Lin Kuo
Journal: International Journal of Aerospace Engineering
 (2016)

---

## Abstract

This paper presents the nonlinear dynamic modeling and control of a tethered satellite system (TSS), and the control strategy is based on the state-dependent Riccati equation (SDRE). The TSS is modeled by a two-piece dumbbell model, which leads to a set of five nonlinear coupled ordinary differential equations. Two sets of equations of motion are proposed, which are based on the first satellite and the mass center of the TSS. There are two reasons to formulate the two sets of equations. One is to facilitate their mutual comparison due to the complex formulations. The other is to provide them for different application situations. Based on the proposed models, the nonlinear dynamic analysis is performed by numerical simulations. Besides, to reduce the convergence time of the librations of the TSS, the SDRE control with a prescribed degree of stability is developed, and the illustrative examples validate the proposed approach.


---

## Body

## 1. Introduction

A tethered satellite system (TSS) usually consists of several spacecraft, a long tether which connects these spacecraft and a tether deployment/retrieval system. In 1895, Tsiolkovsky conceived the concept of a TSS [1]. In 1974, Colombo et al. proposed a “Shuttle-borne Skyhook” for low-orbital altitude research, which marked the advent of tethered satellite systems (TSS) [2]. Since then, there are several potential applications in space technologies, such as scientific experiments, deployment or retrieval of satellites, power generation, measurement of aerodynamic forces.A classical TSS consists of two satellites connected by a tether, and one of the simplest models is a one-piece dumbbell, where the two satellites are assumed to be point masses and the tether is considered to be a straight massless line [3]. He et al. studied the stability of a TSS by utilizing Taylor expansion formula and stability criterion, and a range-rate control algorithm is developed to achieve the stable control of tether’s deploying, keeping, and retrieving [4]. Nakanishi et al. investigated the in-plane periodic solutions for a dumbbell model in elliptic orbits using bifurcation, and their trajectories are projected on the van der Pol planes to predict when the delayed feedback control will need to act to maintain the periodic motions [5]. Kojima et al. conducted an experimental study to investigate the effect of gravity for the libration of tethered satellite system by delayed feedback control [6]. Zhang et al. presented several criteria on the existence of periodic solutions for a TSS in an elliptical orbit, and the uniqueness of periodic solutions for the TSS in a circular orbit is presented on the basis of coincidence degree theory [7]. Hu et al. presented the theoretical and experimental studies of the deployment and retrieval control of a TSS, and the ideas of online optimization and receding horizon control were applied to design a feedback controller for the TSS [8]. Jung et al. presented a dynamic analysis of a TSS with a moving mass based on a two-piece dumbbell model, where the moving mass is conveyed between two satellites along a straight massless tether [9]. The same authors analyzed the nonlinear dynamic behavior of deployment and retrieval of a three-body TSS on a variable-radius orbit, and the system is modeled as a two-piece dumbbell model with six degrees of freedom, which consists of three point masses and two straight massless tethers [10]. Hong and Varatharajoo discussed the development of mathematical models for a flexible TSS in both planar and coplanar states, where the TSS consists of three rigid bodies with two flexible tethers, each connecting two rigid bodies with one located in the center and serving as the mothership [11]. Hong et al. developed the analytic solutions for a TSS subjected to internal tether tension moment and external aerodynamic torque for spin-up and spin-down maneuvers by the approximation of Euler’s equations of motion via Fresnel integrals and sine and cosine integrals [12]. Furthermore, there are some review papers addressing the technologies and applications of a TSS [13–17].For nonlinear optimal control, Pearson presented a linear time and state-dependent approximation to a nonlinear and nonstationary system optimized with respect to a quadratic performance index by treating it as an instantaneously linear stationary system in 1962 [18]. This approach was further investigated by Coultier et al., who extended the LQR approach to propose a so-called state-dependent Riccati equation (SDRE) method for the nonlinear quadratic regulation problems [19]. Zhao and Deng investigated the relationship between the technique by SDRE and Hamilton-Jacobi-Isaacs equations for nonlinear H∞ control design by establishing Lyapunov matrix equations for partial derivatives of the solution of the SDREs and introducing symmetry measure for some related matrices [20]. Nekoo investigated a general case which has control nonlinearities and time-varying weighting matrix Q is solved with three approaches: exact solution, online control update, and power series approximation [21]. Chang and Bentsman extended the continuous time SDRE technique to discrete time under input and state constraints, yielding constrained discrete time SDRE, referred to as CD-SDRE [22]. Babaei and Salamci presented a model reference adaptive control algorithm for a class of nonlinear systems by using a SDRE technique [23]. Nekoo and Geranmehr used a combination of an SDRE controller and SDRE observer to control a class of nonlinear nonaffine control systems [24]. Wang et al. presented a generalized SDRE approach for continuous time nonlinear systems to achieve a mixed nonlinear quadratic regulator and H∞ control performance criteria [25]. Korayem and Nekoo used the finite-time horizon SDRE for controlling a class of nonlinear systems, and the derivations of the SDRE for two classes of systems were presented in the finite-time horizon [26]. Korayem and Nekoo investigated the finite-time optimal and suboptimal controls for time-varying systems with state and control nonlinearities, and a finite-time constraint imposed on the equation leads to a state-dependent differential Riccati equation (SDDRE) [27]. Furthermore, there are some review papers addressing the SDRE method [28–32].This paper investigated a TSS modeled by a two-piece dumbbell model, which leads to five nonlinear coupled ordinary differential equations. There are two sets of equations of motion formulated. One is based on the first satellite and the other is based on the mass center of the TSS. There are two reasons to formulate the two sets of equations. One is to facilitate their mutual comparison due to the complex formulations. The other is to provide them for different application situations. In order to suppress the librations of the system, the SDRE control is applied, which aims to control a nonlinear system and provides the flexibility of controller design. Besides, to increase the convergence rates of the librations, the prescribed degree of stability is incorporated into the SDRE control. This paper is organized as follows. Section2 formulates the two sets of equations of motion of the TSS. Section 3 presents the SDRE control with the prescribed degree of stability. Section 4 analyzes the nonlinear dynamics of the TSS. Section 5 demonstrates the libration suppressions by the SDRE control. Section 6 summarizes some significant conclusions.


## 2. Dynamic Modeling

A tethered satellite system consists of three satellites connected by a tether, where each satellite is considered as a point mass. Besides, the tether is assumed to be massless, and its length is constant. Thus, the system can be simplified as a two-piece dumbbell model shown in Figure1, where the three points A, B, and D represent three satellites, and they are connected to each other by a constant-length massless tether. These satellites are assumed to be point masses and to be expressed as mA, mB, and mD, which refer to the points A, B, and D, respectively. Besides, there is a fixed coordinate system X-Y defined on the center of Earth O. In this paper, one neglects the environmental disturbance torques due to the aerodynamics, the solar radiation, the Earth magnetic filed, and so forth.Figure 1
A tethered satellite system.This section presents two dynamic models to describe the motions of the tethered satellite system. The formulation of the first one is based on satelliteA, and the other is based on the mass center of the system. Since both formulations describe the same system, they must have the same degrees of freedom, which is equal to five in the system. However, both formulations are expressed in terms of different variables. Each formulation can be used to perform the analysis for a given specific prescribed motion. For instance, if a circular motion is imposed on satellite A, the first set of formulations should be adapted. If the same circular motion is imposed on the mass center of the system, the second set of formulations should be selected.


### 2.1. Formulations Based on PointA

Figure2(a) shows the coordinate system defined at point A, and the position of point A is expressed as the orbital radius vector R⃑A with the true anomaly ψA. Besides, another coordinate system fixed on pointA is defined by the two unit vectors e⃑RA and e⃑ψA. The relative position vectors from points A to D and from D to B are, respectively, expressed as L⃑1 and L⃑2, which can show the configuration of the tether. Referring to Figure 2(a), angles ϕ1 and ϕ2, respectively, represent the directions of L⃑1 and L⃑2 based on the unit vector e⃑RA. The two angles are introduced to investigate the libration of the tethered satellite system based on different coordinate systems. If they have an equal value, the system can be treated as a one-piece dumbbell model.Figure 2
Simplified tethered satellite systems based on (a) pointA and (b) point C.
(a) (b)Referring to Figure2(a), the position vectors of the three point masses from the origin of the X-Y coordinate system are given as(1)R⃑A=RAe⃑RA,R⃑B=RA+L1c1+L2c2e⃑RA+L1s1+L2s2e⃑ψA,R⃑D=RA+L1c1e⃑RA+L1s1e⃑ψA,where RA, RB, RD, L1, and L2 represent the length of vectors R⃑A, R⃑B, R⃑D, L⃑1, and L⃑2, respectively; c1 and c2 refer to cos⁡ϕ1 and cos⁡ϕ2, respectively. Note that the three position vectors are expressed in terms of e⃑RA and e⃑ψA.By differentiating (1) with respect to time, the velocities of the three vectors are given as(2)R⃑˙A=R˙Ae⃑RA+RAψ˙Ae⃑ψA,R⃑˙B=R˙A+L˙1c1-c2-L1ψ˙A+ϕ˙1s1-L2ψ˙A+ϕ˙2s2e⃑RA+RAψ˙A+L˙1s1-s2+L1ψ˙A+ϕ˙1c1+L2ψ˙A+ϕ˙2c2e⃑ψA,R⃑˙D=R˙A+L˙1c1-L1ψ˙A+ϕ˙1s1e⃑RA+RAψ˙A+L˙1s1+L1ψ˙A+ϕ˙1c1e⃑ψA.Therefore, the kinetic energy and the potential energy due to the Earth’s gravitational field are, respectively, written as(3)T=12mAR⃑˙A2+12mBR⃑˙B2+12mDR⃑˙D2,(4)V=-GMemARA-1-GMemBRB-1-GMemDRD-1,where G and Me are the universal gravitational constant and the Earth mass, respectively.The nonlinear dynamic equations of motion of the tethered satellite system can be obtained by applying Euler-Lagrange’s equation as(5)ddt∂Lg∂q˙i-∂Lg∂qi=Qi,where qi is the generalized coordinate, Qi is the nonconservative forces corresponding to the generalized coordinate, and Lg is the Lagrangian as(6)Lg=T-Vin which T and V are the kinetic energy and the potential energy, respectively. They are expressed in (3) and (4).Based on the two-piece dumbbell model, the generalized coordinates areRA, ψA, ϕ1, ϕ2, and L1. Besides, the nonconservative forces Qi refer to the driving forces or torques applied to some points and the tether. Substituting (1) to (4) and (6) into (5), one has the nonlinear dynamic equations of motion as follows.Forqi=RA, (5) leads to(7)mA+mB+mDR¨A-mB+mDL1s1+mBL2s2ψ¨A-mB+mDL1s1ϕ¨1-mBL2s2ϕ¨2+mB+mDc1-mBc2L¨1-mA+mB+mDRAψ˙A2-mB+mD2L˙1ψ˙A+ϕ˙1s1+L1ψ˙A+ϕ˙12c1+mB2L˙1ψ˙A+ϕ˙2s2-L2ψ˙A+ϕ˙22c2+GMemARA-2+GMemBRB-3RBx+GMemDRD-3RDx=FA,where FA is the external force applied to point A. One has(8)Ri=Rix2+Riy21/2,fori=A,B,D,(9)RAx,RAy=RA,0,(10)RBx,RBy=RA+L1c1-c2+Lc2,L1s1-s2+Ls2,(11)RCx,RCy=RA+L1c1,L1s1.Note that (9) to (11) are expressed based on the coordinate system, e⃑RA and e⃑ψA.Forqi=ψA, (5) leads to(12)-mB+mDL1s1-mBL2s2R¨A+mA+mB+mDRA2+mB+mDL12+2RAL1c1+mBL22+2L1L2c12+2RAL2c2ψ¨A+mB+mDL12+RAL1c1+mBL1L2c12ϕ¨1+mBL22+L1L2c12+RAL2c2ϕ¨2+mB+mDRAs1+mB-RAs2+Ls12L¨1+2mA+mB+mDRAR˙Aψ˙A+mB+mD-RAL12ψ˙A+ϕ˙1ϕ˙1s1+2L1L˙1ψ˙A+ϕ˙1+2RAL˙1ψ˙A+ϕ˙1c1+2L1R˙Aψ˙Ac1+mB-2L2L˙1ψ˙A+ϕ˙2+L2-L1L˙12ψ˙A+ϕ˙1+ϕ˙2c12-LL22ψ˙A+ϕ˙1+ϕ˙2ϕ˙1-ϕ˙2s12+LL˙1ϕ˙1-ϕ˙2c12+2L2R˙Aψ˙Ac2-2RAL˙1ψ˙A+ϕ˙2c2-RAL22ψ˙A+ϕ˙2ϕ˙2s2=TO,where TO is the external torque applied to point O and c12 and s12 refer to cos⁡ϕ1-ϕ2 and sin⁡ϕ1-ϕ2, respectively.Forqi=ϕ1, (5) leads to(13)-mB+mDL1s1R¨A+mB+mDL12+L1RAc1+mBL1L2c12ψ¨A+mB+mDL12ϕ¨1+mBL1L2c12ϕ¨2+mBL1s12L¨1+mB+mD2L1L˙1ψ˙A+ϕ˙1+RAL1s1ψ˙A2+2L1R˙Aψ˙Ac1+mB-2L1L˙1ψ˙A+θ˙2c12+L1L2ψ˙A+θ˙22s12+GMemBRB-3-RBxL1s1+RByL1c1+GMemDRD-3-RDxL1s1+RDyL1c1=TA,where TA is the external torque applied to point A.Forqi=ϕ2, (5) leads to(14)-mBL2s2R¨A+mBL22+L1L2c12+RAL2c2ψ¨A+mBL1L2c12ϕ¨1+mBL22ϕ¨2+m2L2s12L¨1+mBRAL2s2ψ˙A2-2L2L˙1ψ˙A+ϕ˙2+2L2R˙Aψ˙Ac2+2L2c12L˙1ψ˙A+ϕ˙1-L1L2ψ˙A+ϕ˙12s12+GMemBRB-3-RBxL2s2+RByL2c2=TD,where TD is the external torque applied to point D.Forqi=L1, (5) leads to(15)mB+mDc1-mBc2R¨A+mB+mDRAs1+mB-RAs2+L1s12ψ¨A+mBL1s12ϕ¨1+mBL2s12ϕ¨2+mB+mD+mB1-2c12L¨1+mB+mD2R˙Aψ˙As1-RAψ˙A2c1-L1ψ˙A+ϕ˙12+mB-2R˙Aψ˙As2+RAψ˙A2c2+L2ψ˙A+ϕ˙22+L1ψ˙A+ϕ˙12c12-L2ψ˙A+ϕ˙22c12+2L˙1ϕ˙1-ϕ˙2s12+GMemBRB-3RBxc1-c2+RBys1-s2+GMemDRD-3RDxc1+RDys1=FL,where FL is the external force applied to the tether.Equations (7) and (12) to (15) are the nonlinear dynamic equations of motion for the tethered satellite system based on the coordinate system e⃑RA and e⃑ψA.


### 2.2. Formulations Based on the Mass Center

Figure2(b) shows the coordinate system defined at the mass center, where point C represents the mass center of the tethered satellite system, and the position vector from points O to C is expressed as R⃑C with the true anomaly ψC. The relative position vectors from points C to the three point masses are r⃑A, r⃑B, and r⃑D. Referring to Figure 2(b), another coordinate system fixed on point C is defined by the two unit vectors e⃑RC and e⃑ψC. Referring to Figure 2(b), angles θ1 and θ2, respectively, represent the directions of L⃑1 and L⃑2 based on the unit vector e⃑RC. Since the total length of the tether is a constant value, it can be written as(16)L=L1+L2,where L1 and L2 are the lengths of vectors L⃑1 and L⃑2, respectively.To formulate the dynamic equations of the two-piece dumbbell model, one will apply Euler-Lagrange’s equation, and the position and velocity of each point mass will be expressed in terms of the relative vectorsL⃑1 and L⃑2. Furthermore, they will be written based on the coordinate system (e⃑RC,e⃑ψC).Referring to Figure2(b), one has these two equations as(17)-r⃑A+r⃑D=L⃑1,r⃑B-r⃑D=L⃑2.Besides, since point C is the mass center of the tethered satellite system, one has the equation as(18)ρAr⃑A+ρBr⃑B+ρDr⃑D=0,where the mass ratios are written as(19)ρA=mAM,ρB=mBM,ρD=mDMin which the total mass of the system is given by(20)M=mA+mB+mD.Simultaneously solving (17) and (18) leads to(21)r⃑A=-1-ρAL⃑1-ρBL⃑2,r⃑B=ρAL⃑1+1-ρBL⃑2,r⃑D=ρAL⃑1-ρBL⃑2.Note that (21) expresses the relative vectors r⃑A, r⃑B, and r⃑D in terms of the relative vectors L⃑1 and L⃑2.By taking the time derivatives of (21), the velocities of the relative vectors r⃑A, r⃑B, and r⃑D are given as(22)r⃑˙A=-1-ρAL⃑˙1-ρBL⃑˙2,r⃑˙B=ρAL⃑˙1+1-ρBL⃑˙2,r⃑˙A=ρAL⃑˙1-ρBL⃑˙2.Referring to Figure2(b) again, the relative vectors L⃑1 and L⃑2 can be expressed in the coordinate system (e⃑r,e⃑ψ) as(23)L⃑1=L1c1e⃑RC+L1s1e⃑ψC,L⃑2=L2c2e⃑RC+L2s2e⃑ψC,where c1, c2, s1, and s2 refer to cos⁡θ1, cos⁡θ2, sin⁡θ1, and sin⁡θ2, respectively.By taking the time derivatives of (23), the velocity of the relative vectors L⃑1 and L⃑2 is expressed as(24)L⃑˙1=L˙1c1-L1θ˙1s1-L1ψ˙Cs1e⃑RC+L˙1s1+L1θ˙1c1+L1ψ˙Cc1e⃑ψC,L⃑˙2=L˙2c2-L2θ˙2ss-L2ψ˙Cs2e⃑RC+L˙2s2+L2θ˙2c2+L2ψ˙Cc2e⃑ψC.To obtain the dynamic equations of the system by applying Euler-Lagrange’s equation, the kinetic energy is written as(25)T=12mAR⃑˙C+r⃑˙A2+12mBR⃑˙C+r⃑˙B2+12mDR⃑˙C+r⃑˙D2=12MR⃑˙C2+12mAr⃑˙A2+12mBr⃑˙B2+12mDr⃑˙D2=12MR⃑˙C2+12a1L⃑˙12+12a2L⃑˙22+a3L⃑˙1·L⃑˙2,where(26)a1=mA1-ρA2+mBρA2+mDρA2,a2=mAρB2+mB1-ρB2+mDρB2,a3=mA1-ρAρB+mBρA1-ρB-mDρAρB.And the potential energy is written as(27)V=-GMemAR⃑C+r⃑A-1-GMemBR⃑C+r⃑B-1-GMemDR⃑C+r⃑D-1.Using Euler-Lagrange’s equation shown in (5), one obtains the equations of motion as follows.Forqi=RC, (5) leads to(28)MR¨C-MRCψ˙C2+GMemARA-3RAx+GMemBRB-3RBx+GMemDRD-3RDx=FC,where FC is the external force applied to the center of the system and(29)Ri=Rix2+Riy21/2,fori=A,B,D,RAx,RAy=RC-1-ρAL1c1-ρBL2c2,-1-ρAL1s1-ρBL2s2,RBx,RBy=RC+ρAL1c1+1-ρBL2c2,ρAL1s1+1-ρBL2s2,RCx,RCy=RC+ρAL1c1-ρBL2c2,ρAL1s1-ρBL2s2.Forqi=ψC, (5) leads to(30)MRC2+a1L12+a2L22+2a3L1L2c12ψ¨C+a1L12+a3L1L2c12θ¨1+a2L22+a3L1L2c12θ¨2+a3Ls12L¨1+2MRCR˙Cψ˙C+2a1L1L˙1ψ˙C+θ˙1-2a2L2L˙1ψ˙C+θ˙2+a3L2-L1c12L˙12ψ˙C+θ˙1+θ˙2-a3L1L2s122ψ˙C+θ˙1+θ˙2θ˙1-θ˙2+a3LL˙1c12θ˙1-θ˙2=TO,where TO is the external torque at point O, c12=cos⁡(θ1-θ2), and s12=sin⁡(θ1-θ2).Forqi=θ1, (5) leads to (31)a1L12+a3L1L2c12ψ¨C+a1L12θ¨1+a3L1L2c12θ¨2+a3L1s12L¨1+2a1L1L˙1ψ˙C+θ˙1+a3L1L2c12ψ˙C+θ˙22s12-2a3L1L˙1ψ˙C+θ˙2c12+GMemARA-31-ρAL1s1RAx-c1RAy+GMemBRB-3ρAL1-s1RBx+c1RBy+GMemDRD-3ρAL1-s1RDx+c1RDy=TA,where TA is the external torque applied at point A.Forqi=θ2, (5) leads to(32)a2L22+a3L1L2c12ψ¨C+a3L1L2c12θ¨1+a2L22θ¨2+a3L2s12L¨1-2a2L2L˙1ψ˙C+θ˙2-a3L1L2c12ψ˙C+θ˙12s12+2a3L2L˙1ψ˙C+θ˙1c12+GMemARA-3ρBL2s2RAx-c2RAy+GMemBRB-31-ρBL2-s2RBx+c2RBy+GMemDRD-3ρBL2s2RDx-c2RDy=TD,where TD is the external torque applied at point D.Forqi=L1, (5) leads to(33)a3Ls12ψ¨C+a3L1s12θ¨1+a3L2s12θ¨2+a1+a2-2a3c12L¨1-a1L1ψ˙C+θ˙12+a2L2ψ˙C+θ˙22-a3L2-L1ψ˙C+θ˙1ψ˙C+θ˙2c12+2a3s12L˙1θ˙1-θ˙2+a3L1c12ψ˙C+θ˙1θ˙1-θ˙2+a3L2c12ψ˙C+θ˙2θ˙1-θ˙2+GMemARA-3-1-ρAc1+ρBc2RAx+-1-ρAs1+ρBs2RAy+GMemBRB-3ρAc1-1-ρBc2RBx+ρAs1-1-ρBs2RBy+GMemDRD-3ρAc1+ρBc2RDx+ρAs1+ρBs2RDy=FL,where FL is the external force applied at the tether.Equations (28) and (30) to (33) are the nonlinear dynamic equations of motion for the tethered satellite system based on the coordinate system (e⃑RC,e⃑ψC). Note that the term related to R¨C only appears in (28) but not in (30) to (33). This property is not shown in the equations of motion based on point A.


## 2.1. Formulations Based on PointA

Figure2(a) shows the coordinate system defined at point A, and the position of point A is expressed as the orbital radius vector R⃑A with the true anomaly ψA. Besides, another coordinate system fixed on pointA is defined by the two unit vectors e⃑RA and e⃑ψA. The relative position vectors from points A to D and from D to B are, respectively, expressed as L⃑1 and L⃑2, which can show the configuration of the tether. Referring to Figure 2(a), angles ϕ1 and ϕ2, respectively, represent the directions of L⃑1 and L⃑2 based on the unit vector e⃑RA. The two angles are introduced to investigate the libration of the tethered satellite system based on different coordinate systems. If they have an equal value, the system can be treated as a one-piece dumbbell model.Figure 2
Simplified tethered satellite systems based on (a) pointA and (b) point C.
(a) (b)Referring to Figure2(a), the position vectors of the three point masses from the origin of the X-Y coordinate system are given as(1)R⃑A=RAe⃑RA,R⃑B=RA+L1c1+L2c2e⃑RA+L1s1+L2s2e⃑ψA,R⃑D=RA+L1c1e⃑RA+L1s1e⃑ψA,where RA, RB, RD, L1, and L2 represent the length of vectors R⃑A, R⃑B, R⃑D, L⃑1, and L⃑2, respectively; c1 and c2 refer to cos⁡ϕ1 and cos⁡ϕ2, respectively. Note that the three position vectors are expressed in terms of e⃑RA and e⃑ψA.By differentiating (1) with respect to time, the velocities of the three vectors are given as(2)R⃑˙A=R˙Ae⃑RA+RAψ˙Ae⃑ψA,R⃑˙B=R˙A+L˙1c1-c2-L1ψ˙A+ϕ˙1s1-L2ψ˙A+ϕ˙2s2e⃑RA+RAψ˙A+L˙1s1-s2+L1ψ˙A+ϕ˙1c1+L2ψ˙A+ϕ˙2c2e⃑ψA,R⃑˙D=R˙A+L˙1c1-L1ψ˙A+ϕ˙1s1e⃑RA+RAψ˙A+L˙1s1+L1ψ˙A+ϕ˙1c1e⃑ψA.Therefore, the kinetic energy and the potential energy due to the Earth’s gravitational field are, respectively, written as(3)T=12mAR⃑˙A2+12mBR⃑˙B2+12mDR⃑˙D2,(4)V=-GMemARA-1-GMemBRB-1-GMemDRD-1,where G and Me are the universal gravitational constant and the Earth mass, respectively.The nonlinear dynamic equations of motion of the tethered satellite system can be obtained by applying Euler-Lagrange’s equation as(5)ddt∂Lg∂q˙i-∂Lg∂qi=Qi,where qi is the generalized coordinate, Qi is the nonconservative forces corresponding to the generalized coordinate, and Lg is the Lagrangian as(6)Lg=T-Vin which T and V are the kinetic energy and the potential energy, respectively. They are expressed in (3) and (4).Based on the two-piece dumbbell model, the generalized coordinates areRA, ψA, ϕ1, ϕ2, and L1. Besides, the nonconservative forces Qi refer to the driving forces or torques applied to some points and the tether. Substituting (1) to (4) and (6) into (5), one has the nonlinear dynamic equations of motion as follows.Forqi=RA, (5) leads to(7)mA+mB+mDR¨A-mB+mDL1s1+mBL2s2ψ¨A-mB+mDL1s1ϕ¨1-mBL2s2ϕ¨2+mB+mDc1-mBc2L¨1-mA+mB+mDRAψ˙A2-mB+mD2L˙1ψ˙A+ϕ˙1s1+L1ψ˙A+ϕ˙12c1+mB2L˙1ψ˙A+ϕ˙2s2-L2ψ˙A+ϕ˙22c2+GMemARA-2+GMemBRB-3RBx+GMemDRD-3RDx=FA,where FA is the external force applied to point A. One has(8)Ri=Rix2+Riy21/2,fori=A,B,D,(9)RAx,RAy=RA,0,(10)RBx,RBy=RA+L1c1-c2+Lc2,L1s1-s2+Ls2,(11)RCx,RCy=RA+L1c1,L1s1.Note that (9) to (11) are expressed based on the coordinate system, e⃑RA and e⃑ψA.Forqi=ψA, (5) leads to(12)-mB+mDL1s1-mBL2s2R¨A+mA+mB+mDRA2+mB+mDL12+2RAL1c1+mBL22+2L1L2c12+2RAL2c2ψ¨A+mB+mDL12+RAL1c1+mBL1L2c12ϕ¨1+mBL22+L1L2c12+RAL2c2ϕ¨2+mB+mDRAs1+mB-RAs2+Ls12L¨1+2mA+mB+mDRAR˙Aψ˙A+mB+mD-RAL12ψ˙A+ϕ˙1ϕ˙1s1+2L1L˙1ψ˙A+ϕ˙1+2RAL˙1ψ˙A+ϕ˙1c1+2L1R˙Aψ˙Ac1+mB-2L2L˙1ψ˙A+ϕ˙2+L2-L1L˙12ψ˙A+ϕ˙1+ϕ˙2c12-LL22ψ˙A+ϕ˙1+ϕ˙2ϕ˙1-ϕ˙2s12+LL˙1ϕ˙1-ϕ˙2c12+2L2R˙Aψ˙Ac2-2RAL˙1ψ˙A+ϕ˙2c2-RAL22ψ˙A+ϕ˙2ϕ˙2s2=TO,where TO is the external torque applied to point O and c12 and s12 refer to cos⁡ϕ1-ϕ2 and sin⁡ϕ1-ϕ2, respectively.Forqi=ϕ1, (5) leads to(13)-mB+mDL1s1R¨A+mB+mDL12+L1RAc1+mBL1L2c12ψ¨A+mB+mDL12ϕ¨1+mBL1L2c12ϕ¨2+mBL1s12L¨1+mB+mD2L1L˙1ψ˙A+ϕ˙1+RAL1s1ψ˙A2+2L1R˙Aψ˙Ac1+mB-2L1L˙1ψ˙A+θ˙2c12+L1L2ψ˙A+θ˙22s12+GMemBRB-3-RBxL1s1+RByL1c1+GMemDRD-3-RDxL1s1+RDyL1c1=TA,where TA is the external torque applied to point A.Forqi=ϕ2, (5) leads to(14)-mBL2s2R¨A+mBL22+L1L2c12+RAL2c2ψ¨A+mBL1L2c12ϕ¨1+mBL22ϕ¨2+m2L2s12L¨1+mBRAL2s2ψ˙A2-2L2L˙1ψ˙A+ϕ˙2+2L2R˙Aψ˙Ac2+2L2c12L˙1ψ˙A+ϕ˙1-L1L2ψ˙A+ϕ˙12s12+GMemBRB-3-RBxL2s2+RByL2c2=TD,where TD is the external torque applied to point D.Forqi=L1, (5) leads to(15)mB+mDc1-mBc2R¨A+mB+mDRAs1+mB-RAs2+L1s12ψ¨A+mBL1s12ϕ¨1+mBL2s12ϕ¨2+mB+mD+mB1-2c12L¨1+mB+mD2R˙Aψ˙As1-RAψ˙A2c1-L1ψ˙A+ϕ˙12+mB-2R˙Aψ˙As2+RAψ˙A2c2+L2ψ˙A+ϕ˙22+L1ψ˙A+ϕ˙12c12-L2ψ˙A+ϕ˙22c12+2L˙1ϕ˙1-ϕ˙2s12+GMemBRB-3RBxc1-c2+RBys1-s2+GMemDRD-3RDxc1+RDys1=FL,where FL is the external force applied to the tether.Equations (7) and (12) to (15) are the nonlinear dynamic equations of motion for the tethered satellite system based on the coordinate system e⃑RA and e⃑ψA.


## 2.2. Formulations Based on the Mass Center

Figure2(b) shows the coordinate system defined at the mass center, where point C represents the mass center of the tethered satellite system, and the position vector from points O to C is expressed as R⃑C with the true anomaly ψC. The relative position vectors from points C to the three point masses are r⃑A, r⃑B, and r⃑D. Referring to Figure 2(b), another coordinate system fixed on point C is defined by the two unit vectors e⃑RC and e⃑ψC. Referring to Figure 2(b), angles θ1 and θ2, respectively, represent the directions of L⃑1 and L⃑2 based on the unit vector e⃑RC. Since the total length of the tether is a constant value, it can be written as(16)L=L1+L2,where L1 and L2 are the lengths of vectors L⃑1 and L⃑2, respectively.To formulate the dynamic equations of the two-piece dumbbell model, one will apply Euler-Lagrange’s equation, and the position and velocity of each point mass will be expressed in terms of the relative vectorsL⃑1 and L⃑2. Furthermore, they will be written based on the coordinate system (e⃑RC,e⃑ψC).Referring to Figure2(b), one has these two equations as(17)-r⃑A+r⃑D=L⃑1,r⃑B-r⃑D=L⃑2.Besides, since point C is the mass center of the tethered satellite system, one has the equation as(18)ρAr⃑A+ρBr⃑B+ρDr⃑D=0,where the mass ratios are written as(19)ρA=mAM,ρB=mBM,ρD=mDMin which the total mass of the system is given by(20)M=mA+mB+mD.Simultaneously solving (17) and (18) leads to(21)r⃑A=-1-ρAL⃑1-ρBL⃑2,r⃑B=ρAL⃑1+1-ρBL⃑2,r⃑D=ρAL⃑1-ρBL⃑2.Note that (21) expresses the relative vectors r⃑A, r⃑B, and r⃑D in terms of the relative vectors L⃑1 and L⃑2.By taking the time derivatives of (21), the velocities of the relative vectors r⃑A, r⃑B, and r⃑D are given as(22)r⃑˙A=-1-ρAL⃑˙1-ρBL⃑˙2,r⃑˙B=ρAL⃑˙1+1-ρBL⃑˙2,r⃑˙A=ρAL⃑˙1-ρBL⃑˙2.Referring to Figure2(b) again, the relative vectors L⃑1 and L⃑2 can be expressed in the coordinate system (e⃑r,e⃑ψ) as(23)L⃑1=L1c1e⃑RC+L1s1e⃑ψC,L⃑2=L2c2e⃑RC+L2s2e⃑ψC,where c1, c2, s1, and s2 refer to cos⁡θ1, cos⁡θ2, sin⁡θ1, and sin⁡θ2, respectively.By taking the time derivatives of (23), the velocity of the relative vectors L⃑1 and L⃑2 is expressed as(24)L⃑˙1=L˙1c1-L1θ˙1s1-L1ψ˙Cs1e⃑RC+L˙1s1+L1θ˙1c1+L1ψ˙Cc1e⃑ψC,L⃑˙2=L˙2c2-L2θ˙2ss-L2ψ˙Cs2e⃑RC+L˙2s2+L2θ˙2c2+L2ψ˙Cc2e⃑ψC.To obtain the dynamic equations of the system by applying Euler-Lagrange’s equation, the kinetic energy is written as(25)T=12mAR⃑˙C+r⃑˙A2+12mBR⃑˙C+r⃑˙B2+12mDR⃑˙C+r⃑˙D2=12MR⃑˙C2+12mAr⃑˙A2+12mBr⃑˙B2+12mDr⃑˙D2=12MR⃑˙C2+12a1L⃑˙12+12a2L⃑˙22+a3L⃑˙1·L⃑˙2,where(26)a1=mA1-ρA2+mBρA2+mDρA2,a2=mAρB2+mB1-ρB2+mDρB2,a3=mA1-ρAρB+mBρA1-ρB-mDρAρB.And the potential energy is written as(27)V=-GMemAR⃑C+r⃑A-1-GMemBR⃑C+r⃑B-1-GMemDR⃑C+r⃑D-1.Using Euler-Lagrange’s equation shown in (5), one obtains the equations of motion as follows.Forqi=RC, (5) leads to(28)MR¨C-MRCψ˙C2+GMemARA-3RAx+GMemBRB-3RBx+GMemDRD-3RDx=FC,where FC is the external force applied to the center of the system and(29)Ri=Rix2+Riy21/2,fori=A,B,D,RAx,RAy=RC-1-ρAL1c1-ρBL2c2,-1-ρAL1s1-ρBL2s2,RBx,RBy=RC+ρAL1c1+1-ρBL2c2,ρAL1s1+1-ρBL2s2,RCx,RCy=RC+ρAL1c1-ρBL2c2,ρAL1s1-ρBL2s2.Forqi=ψC, (5) leads to(30)MRC2+a1L12+a2L22+2a3L1L2c12ψ¨C+a1L12+a3L1L2c12θ¨1+a2L22+a3L1L2c12θ¨2+a3Ls12L¨1+2MRCR˙Cψ˙C+2a1L1L˙1ψ˙C+θ˙1-2a2L2L˙1ψ˙C+θ˙2+a3L2-L1c12L˙12ψ˙C+θ˙1+θ˙2-a3L1L2s122ψ˙C+θ˙1+θ˙2θ˙1-θ˙2+a3LL˙1c12θ˙1-θ˙2=TO,where TO is the external torque at point O, c12=cos⁡(θ1-θ2), and s12=sin⁡(θ1-θ2).Forqi=θ1, (5) leads to (31)a1L12+a3L1L2c12ψ¨C+a1L12θ¨1+a3L1L2c12θ¨2+a3L1s12L¨1+2a1L1L˙1ψ˙C+θ˙1+a3L1L2c12ψ˙C+θ˙22s12-2a3L1L˙1ψ˙C+θ˙2c12+GMemARA-31-ρAL1s1RAx-c1RAy+GMemBRB-3ρAL1-s1RBx+c1RBy+GMemDRD-3ρAL1-s1RDx+c1RDy=TA,where TA is the external torque applied at point A.Forqi=θ2, (5) leads to(32)a2L22+a3L1L2c12ψ¨C+a3L1L2c12θ¨1+a2L22θ¨2+a3L2s12L¨1-2a2L2L˙1ψ˙C+θ˙2-a3L1L2c12ψ˙C+θ˙12s12+2a3L2L˙1ψ˙C+θ˙1c12+GMemARA-3ρBL2s2RAx-c2RAy+GMemBRB-31-ρBL2-s2RBx+c2RBy+GMemDRD-3ρBL2s2RDx-c2RDy=TD,where TD is the external torque applied at point D.Forqi=L1, (5) leads to(33)a3Ls12ψ¨C+a3L1s12θ¨1+a3L2s12θ¨2+a1+a2-2a3c12L¨1-a1L1ψ˙C+θ˙12+a2L2ψ˙C+θ˙22-a3L2-L1ψ˙C+θ˙1ψ˙C+θ˙2c12+2a3s12L˙1θ˙1-θ˙2+a3L1c12ψ˙C+θ˙1θ˙1-θ˙2+a3L2c12ψ˙C+θ˙2θ˙1-θ˙2+GMemARA-3-1-ρAc1+ρBc2RAx+-1-ρAs1+ρBs2RAy+GMemBRB-3ρAc1-1-ρBc2RBx+ρAs1-1-ρBs2RBy+GMemDRD-3ρAc1+ρBc2RDx+ρAs1+ρBs2RDy=FL,where FL is the external force applied at the tether.Equations (28) and (30) to (33) are the nonlinear dynamic equations of motion for the tethered satellite system based on the coordinate system (e⃑RC,e⃑ψC). Note that the term related to R¨C only appears in (28) but not in (30) to (33). This property is not shown in the equations of motion based on point A.


## 3. SDRE Control with a Specified Degree of Stability

A nonlinear and observable system is written as(34)x˙t=fx+gxut,x0=x0,where x(t) and u(t) are the state and input vectors, respectively.To perform a regulation problem, an infinite-time performance index is written as(35)Jx0,ut=12∫0∞xTxQxxt+uTtRxutdt,where Q(t)≥0 and R(t)>0 are state-dependent weighting matrices.To apply the LQR, (34) can be factorized in a linear-like structure as(36)x˙t=Axxt+Bxut,where f(x)=A(x)x(t) and g(x)=B(x). There is an infinite number of forms to factorizef(x) as A(x)x(t), so this property increases the flexibility of the SDRE.Solving (35) and (36) by following the formulation of the LQR leads to(37)ut=-Kxxt=-R-1xBTxPxxt,where P(x) is the symmetric and positive-definite solution of the algebraic state-dependent Riccati equation as(38)PxAx+ATxPx-PxBxR-1xBTxPx+Qx=0.To ensure the closed-loop system with a specified degree of stability, the performance index is rewritten as(39)Jx0,ut=12∫0∞e2αtxTxQxxt+uTtRxutdt,where α is a positive parameter. Note that the pair A+αI,B is completely stabilizable, where I is an identity matrix.Define new state variables and controls as(40)x^t=eαtxt,u^t=eαtut.Using (40), the system equations shown in (36) can be modified as(41)x^˙t=Ax+αIx^t+Bu^t.Note that the initial condition is the same as x^(0)=x(0)=x0. Also, the performance index is transformed as(42)Jx0,u^t=12∫0∞x^TxQxx^t+u^TtRxu^tdt.Solving (41) and (42) by following the formulation of the LQR leads to(43)u^t=-K-xx^t=-R-1xBTxP-xx^t,where P-(x) is the symmetric and positive-definite solution of the algebraic state-dependent Riccati equation as(44)P-xAx+αI+Ax+αITP-x-P-xBxR-1xBTxP-x+Qx=0.Substituting (43) into (41) leads to the optimal closed-loop system as(45)x^˙t=Ax+αI-BxR-1xBTxP-xx^t.Roughly speaking, the system seemingly has eigenvalues with real parts less than -α, so the system is treated with a degree of stability of at least α.Based on (43), the optimal control can be transformed to(46)ut=e-αtu^t=-K-xx=-R-1xBTxP-xxt.Note that the original control (see (37)) has the same structure as the optimal control (see (46)) of the modified system.


## 4. Nonlinear Dynamic Analysis

This section presents the nonlinear dynamic analysis without applying any control theory.


### 4.1. Free Libration

To demonstrate the consistence of the proposed formulations, the numerical simulations of free libration are performed. The initial configuration is that three satellites and the mass center of the Earth are collinear, and satelliteD is at the middle point between satellites A and B. Figure 3 shows the initial configuration, where the distance from the Earth center to satellite A is 42,522 km and the length of the tether is 10,000 km. The three satellites have the same initial orbital rate ψ˙A=7.2923×10-5 rad/s. Hence, the initial velocities of these satellites can be determined and are expressed as v⃑A, v⃑B, and v⃑D shown in Figure 3. The masses of satellites are identical as 10,000 kg, and the mass of satellite is 100 kg. Based on the aforementioned initial state, the initial conditions are given as RA = 42,522 km, ψA=0, ϕ1=0, ϕ2=0, L1 = 5,000 km, R˙A=0, ψ˙A=7.2923×10-5 rad/s, ϕ˙1=0, ϕ˙2=0, and L˙1=0 for the formulations based on point A. The initial conditions are RC = 47,522 km, ψC=0, θ1=0, θ2=0, L1=5000 km, R˙C=0, ψ˙C=7.2923×10-5 rad/s, θ˙1=0, θ˙2=0, and L˙1=0 for the formulations based on pointC.Figure 3
Initial configuration for free libration.Since both formulations represent the same system and the same initial state is applied, the numerical simulations must be identical. Figure4 shows distances RA and RC as functions to time. The results show that the distances periodically change and one period is around 2.389×104 seconds. Besides, the mean values of RA and RC are 7.7083×104 and 8.1791×104 km, respectively. Also, their amplitude is 3.4561×104 and 3.4269×104 km, respectively. Figure 5 shows the true anomalies ψA and ψC as functions of time. The results show that they do not linearly change with time. Figure 6 shows the libration angles ϕ1 and ϕ2 and the tether length L1 as functions of time. Since the total length of the tether is 10,000 km, the tether length L1 must be less than the total length. This figure shows that length L1 changes between 0 and 10,000 km. Besides, while length L1 reaches maximum value 10,000 km, the libration angle ϕ2 will have the change of 360°. Similarly, while length L1 reaches the minimum value 0, the libration angle ϕ1 also has the change of 360°. Figure 7 shows the libration angles θ1 and θ2 and tether length L1 as functions of time, and this figure is similar to Figure 6.Figure 4
The distances from Earth center as functions of time for free libration.Figure 5
True anomaly as functions of time for free libration.Figure 6
Libration angles and tether lengthL1 by the formulations based on point A.Figure 7
Libration angles and tether lengthL1 by the formulations based on point C.


### 4.2. Free Librations with a Prescribed FunctionL1(t)

One considers that satelliteD is a climber, which will move along the tether based on a prescribed function. The climber is initially located at L10=0.1 km. It ascends while accelerating until reaching L11=0.1L. It cruises at a constant velocity vc until L12=L1f-0.1L, where L1f is the final location of the climber. During the acceleration and deceleration periods, the velocity varies as a sinusoidal function. Hence, the motion of the climber is described by the following piecewise function as(47)L1t=L10+vc2t-t1πsin⁡πtt1,for0≤t≤t1L10+vct-t12,fort1<t≤t2L10+vc2t+t2-t1+tf-t2πsin⁡πt-t2tf-t2,fort2<t≤tfL1f,fort>tf,where(48)t1=2L11-L10vc,t2=L12+L11-2L10vc,tf=2L1f-L10+L11-L12vc.In this section, the initial conditions are the same as those in Section4.1, except for the fact that L1=5000 km and L˙1=0 are removed. They are replaced by the prescribed function shown in (47) to (48), where the constant velocity vc and the final position L1f are specified as 0.5 km/s and 4,628 km, respectively. Since both formulations based on points A and C represent the same system, and the same initial state is considered, the numerical simulations by both formulations must be identical. Figure 8 shows distances RA and RC as functions to time. The results show that the distances periodically change, and one period is around 2.372×104 seconds. Besides, the mean values of RA and RC are 7.6706×104 and 8.1393×104 km, respectively. Also, their amplitude is 3.4184×104 and 3.3896×104 km, respectively. Figure 9 shows the true anomalies ψA and ψC as functions of time. The results show that they do not linearly change with time. To compare the results with those in Section 4.1, all the values are similar, so the prescribed function does not affect RA and RC very much. Figure 10 shows libration angles ϕ1 and ϕ2 and tether length L1 as functions of time. The results show that this prescribed function L1(t) yields the librations angles varying around between ±50°. Figure 11 shows libration angles θ1 and θ2 as well as tether length L1 as functions of time, and the results are similar to those shown in Figure 10.Figure 8
The distances from Earth center as function of time for a prescribed functionL1(t).Figure 9
True anomaly as function of time for a prescribed functionL1(t).Figure 10
Libration angles and tether lengthL1 by the formulations based on point A for a prescribed function L1(t).Figure 11
Libration angles and tether lengthL1 by the formulations based on point C with a prescribed function L1(t).


### 4.3. Free Libration with the Prescribed FunctionL1(t) on a Circular Orbit with a Constant Angular Rate

In this subsection, one considers that the tethered satellite system moves on a circular orbit, so the radius and the true anomaly become known parameters, and the system is governed by two ordinary differential equations. The tether length functionL1(t) is the same as (47).For the formulations based on pointA, RA = 42,522 km, ψ˙A=7.2923×10-5, and R˙A=R¨A=ψA(0)=ψ¨A=0. The initial conditions are set as zeros. Figure 12 shows the libration angles ϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit. The results show that there is a transient response during the first 1.5×105 seconds, and the steady state shows that the libration angles oscillate around zero degrees with the amplitude of 0.01 degrees. For the formulations based on pointC, RC = 42,522 km, ψ˙C=7.2923×10-5, and R˙C=R¨C=ψC(0)=ψ¨C=0. The initial conditions are also set as zeros. Note that the initial state is different from the former case, because the circular orbit is assigned to different points. Figure 13 shows the libration angles ϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit. The results show that there is a transient response during the first 1.5×105 seconds, and the steady state shows that the libration angles oscillate around zero degrees with the amplitude of 0.01 degrees.Figure 12
Libration anglesϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit with zero initial conditions.Figure 13
Libration anglesθ1(t) and θ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit with zero initial conditions.The following simulations are similar to the former ones, except for the fact that the initial conditions are replaced byϕ1=ϕ2=θ1=θ2=0.1 degrees. Based on the formulations based on pointA, Figure 14 shows the libration angles ϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit. The results show that there is a transient response during the first 1.5×105 seconds, and the steady state shows that the libration angles oscillate around zero degrees with the amplitude of 0.1 degrees. Based on the formulations based on pointC, Figure 15 shows the libration angles ϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit. The results show that there is a transient response during the first 1.5×105 seconds, and the steady state shows that the libration angles oscillate around zero degrees with the amplitude of 0.1 degrees.Figure 14
Libration anglesϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit.Figure 15
Libration anglesθ1(t) and θ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit.


## 4.1. Free Libration

To demonstrate the consistence of the proposed formulations, the numerical simulations of free libration are performed. The initial configuration is that three satellites and the mass center of the Earth are collinear, and satelliteD is at the middle point between satellites A and B. Figure 3 shows the initial configuration, where the distance from the Earth center to satellite A is 42,522 km and the length of the tether is 10,000 km. The three satellites have the same initial orbital rate ψ˙A=7.2923×10-5 rad/s. Hence, the initial velocities of these satellites can be determined and are expressed as v⃑A, v⃑B, and v⃑D shown in Figure 3. The masses of satellites are identical as 10,000 kg, and the mass of satellite is 100 kg. Based on the aforementioned initial state, the initial conditions are given as RA = 42,522 km, ψA=0, ϕ1=0, ϕ2=0, L1 = 5,000 km, R˙A=0, ψ˙A=7.2923×10-5 rad/s, ϕ˙1=0, ϕ˙2=0, and L˙1=0 for the formulations based on point A. The initial conditions are RC = 47,522 km, ψC=0, θ1=0, θ2=0, L1=5000 km, R˙C=0, ψ˙C=7.2923×10-5 rad/s, θ˙1=0, θ˙2=0, and L˙1=0 for the formulations based on pointC.Figure 3
Initial configuration for free libration.Since both formulations represent the same system and the same initial state is applied, the numerical simulations must be identical. Figure4 shows distances RA and RC as functions to time. The results show that the distances periodically change and one period is around 2.389×104 seconds. Besides, the mean values of RA and RC are 7.7083×104 and 8.1791×104 km, respectively. Also, their amplitude is 3.4561×104 and 3.4269×104 km, respectively. Figure 5 shows the true anomalies ψA and ψC as functions of time. The results show that they do not linearly change with time. Figure 6 shows the libration angles ϕ1 and ϕ2 and the tether length L1 as functions of time. Since the total length of the tether is 10,000 km, the tether length L1 must be less than the total length. This figure shows that length L1 changes between 0 and 10,000 km. Besides, while length L1 reaches maximum value 10,000 km, the libration angle ϕ2 will have the change of 360°. Similarly, while length L1 reaches the minimum value 0, the libration angle ϕ1 also has the change of 360°. Figure 7 shows the libration angles θ1 and θ2 and tether length L1 as functions of time, and this figure is similar to Figure 6.Figure 4
The distances from Earth center as functions of time for free libration.Figure 5
True anomaly as functions of time for free libration.Figure 6
Libration angles and tether lengthL1 by the formulations based on point A.Figure 7
Libration angles and tether lengthL1 by the formulations based on point C.


## 4.2. Free Librations with a Prescribed FunctionL1(t)

One considers that satelliteD is a climber, which will move along the tether based on a prescribed function. The climber is initially located at L10=0.1 km. It ascends while accelerating until reaching L11=0.1L. It cruises at a constant velocity vc until L12=L1f-0.1L, where L1f is the final location of the climber. During the acceleration and deceleration periods, the velocity varies as a sinusoidal function. Hence, the motion of the climber is described by the following piecewise function as(47)L1t=L10+vc2t-t1πsin⁡πtt1,for0≤t≤t1L10+vct-t12,fort1<t≤t2L10+vc2t+t2-t1+tf-t2πsin⁡πt-t2tf-t2,fort2<t≤tfL1f,fort>tf,where(48)t1=2L11-L10vc,t2=L12+L11-2L10vc,tf=2L1f-L10+L11-L12vc.In this section, the initial conditions are the same as those in Section4.1, except for the fact that L1=5000 km and L˙1=0 are removed. They are replaced by the prescribed function shown in (47) to (48), where the constant velocity vc and the final position L1f are specified as 0.5 km/s and 4,628 km, respectively. Since both formulations based on points A and C represent the same system, and the same initial state is considered, the numerical simulations by both formulations must be identical. Figure 8 shows distances RA and RC as functions to time. The results show that the distances periodically change, and one period is around 2.372×104 seconds. Besides, the mean values of RA and RC are 7.6706×104 and 8.1393×104 km, respectively. Also, their amplitude is 3.4184×104 and 3.3896×104 km, respectively. Figure 9 shows the true anomalies ψA and ψC as functions of time. The results show that they do not linearly change with time. To compare the results with those in Section 4.1, all the values are similar, so the prescribed function does not affect RA and RC very much. Figure 10 shows libration angles ϕ1 and ϕ2 and tether length L1 as functions of time. The results show that this prescribed function L1(t) yields the librations angles varying around between ±50°. Figure 11 shows libration angles θ1 and θ2 as well as tether length L1 as functions of time, and the results are similar to those shown in Figure 10.Figure 8
The distances from Earth center as function of time for a prescribed functionL1(t).Figure 9
True anomaly as function of time for a prescribed functionL1(t).Figure 10
Libration angles and tether lengthL1 by the formulations based on point A for a prescribed function L1(t).Figure 11
Libration angles and tether lengthL1 by the formulations based on point C with a prescribed function L1(t).


## 4.3. Free Libration with the Prescribed FunctionL1(t) on a Circular Orbit with a Constant Angular Rate

In this subsection, one considers that the tethered satellite system moves on a circular orbit, so the radius and the true anomaly become known parameters, and the system is governed by two ordinary differential equations. The tether length functionL1(t) is the same as (47).For the formulations based on pointA, RA = 42,522 km, ψ˙A=7.2923×10-5, and R˙A=R¨A=ψA(0)=ψ¨A=0. The initial conditions are set as zeros. Figure 12 shows the libration angles ϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit. The results show that there is a transient response during the first 1.5×105 seconds, and the steady state shows that the libration angles oscillate around zero degrees with the amplitude of 0.01 degrees. For the formulations based on pointC, RC = 42,522 km, ψ˙C=7.2923×10-5, and R˙C=R¨C=ψC(0)=ψ¨C=0. The initial conditions are also set as zeros. Note that the initial state is different from the former case, because the circular orbit is assigned to different points. Figure 13 shows the libration angles ϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit. The results show that there is a transient response during the first 1.5×105 seconds, and the steady state shows that the libration angles oscillate around zero degrees with the amplitude of 0.01 degrees.Figure 12
Libration anglesϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit with zero initial conditions.Figure 13
Libration anglesθ1(t) and θ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit with zero initial conditions.The following simulations are similar to the former ones, except for the fact that the initial conditions are replaced byϕ1=ϕ2=θ1=θ2=0.1 degrees. Based on the formulations based on pointA, Figure 14 shows the libration angles ϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit. The results show that there is a transient response during the first 1.5×105 seconds, and the steady state shows that the libration angles oscillate around zero degrees with the amplitude of 0.1 degrees. Based on the formulations based on pointC, Figure 15 shows the libration angles ϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit. The results show that there is a transient response during the first 1.5×105 seconds, and the steady state shows that the libration angles oscillate around zero degrees with the amplitude of 0.1 degrees.Figure 14
Libration anglesϕ1(t) and ϕ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit.Figure 15
Libration anglesθ1(t) and θ2(t) as functions of time based on the prescribed function L1(t) on a circular orbit.


## 5. Libration Suppression

This section demonstrates the libration suppression by applying the SDRE control with a prescribed degree of stability, which is presented in Section3.


### 5.1. Libration Control for the Equations of Motion Based on PointA

One considers that the tethered satellite system moves on a circular orbit with a constant rotating rate and satelliteD is fixedly located at the middle point of the tether. Thus, there are only two ordinary differential equations taken into account, and each equation is considered with a control torque, which is obtained based on the SDRE control method.Based on the formulations of the SDRE control, the two equations of motion are factorized as(49)Aϕϕ¨+Bϕϕ˙+Cϕϕ=Du,where A(ϕ), B(ϕ), and C(ϕ) are 2×2 matrices, D is a 2×2 identity matrix, ϕ=ϕ1ϕ2T, and u=TATDT. Besides, the components of the matrices and the vector are written as(50)A11ϕ=mB+mDL12,A12ϕ=A21ϕ=mBL1L2c12,A22ϕ=mBL22,B11ϕ=2mB+mDL1L˙1,B12ϕ=-2mDL1L˙1c12,B21ϕ=2mDL2L˙1c12,B22ϕ=-2mBL2L˙1,C11ϕ=μmBRB-3-RAL1s1+L1L2s1c2+μmDRD-3-RAL1s1+mB+mD2L1L˙1ψ˙A+RAL1ψ˙A2-L1R¨As1+mBL1L2ψ˙A+ϕ˙22s1c2+L1L¨1s1c2ϕ1-1,C12ϕ=μmBRB-3-L1L2c1s2+mB+mD2L1R˙Aψ˙Ac1+L12+RAL1c1ψ¨A+mB-2L1L˙1ψ˙Ac12-L1L2ψ˙A+ϕ˙22c1s2+L1L2ψ¨Ac12-L1L¨1c1s2ϕ2-1,C21ϕ=μmBRB-3L1L2s1c2+mB-2L2L˙1ψ˙Ac12+2L2R˙Aψ˙Ac2-L1L2ψ˙A+ϕ˙12s1c2+L22+L1L2c12+RAL2c2ψ¨A+L2L¨1s1c2ϕ1-1,C22ϕ=μmBRB-3RAL1s2-L12c1s2+mB-2L2L˙1ψ˙A+RAL2ψ˙A2s2+L1L2ψ˙A+ϕ˙12c1s2-L2R¨As2-L2L¨1c1s2ϕ2-1.The simulations are performed based on the following parameters:RA = 42,522 km, ψ˙A=7.2923×10-5 rad/s, and R˙A=R¨A=ψA(0)=ψ¨A=0. Besides, the initial conditions are given as ϕ1=1°, ϕ2=-1°, and ϕ˙1=ϕ˙2=0. For simplicity, the weighting matrices in the SDRE control are selected as Q=diag⁡(1,1,0,0) and R=diag⁡(1,1). In order to demonstrate the performance of the prescribed degree of stability in the SDRE control, there are two cases tested as follows.Case 1 (α=0).
Figure16 shows the libration angles as functions of time, and the results show that the SDRE control can suppress the librations around within 10 orbits. Although the time responses seemingly have high-frequency oscillations in the transient states (see Figure 16(a)), the oscillation period is around 1,930 seconds (see Figure 16(b)). Figure 17 shows the control torques as function of time, and the results show that the maximum control torque is around 0.16 MN-m.Figure 16
Libration angles from the formulations based on pointA by applying the SDRE control with α=0 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Figure 17
Control torques from the formulations based on pointA by applying the SDRE control with α=0 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Case 2 (α=10-6).
Seemingly, the convergence rate forα=0 is little slow, so this test case is based on α=10-6. Figure 18 shows the libration angles as functions of time, and the results show that the SDRE control can suppress the librations around within 8 orbits, which indicates that the proposed SDRE control with a prescribed degree of stability can effectively reduce the converge time. This figure also shows a high-frequency oscillation in the transient state, and the oscillation period is the same as 1,930 seconds. Figure 19 shows the control torques as function of time, and the results show that the maximum control torque is around 0.28 MN-m, which is larger than those presented in previous case. This is the penalty to reduce the convergence time.Figure 18
Libration angles from the formulations based on pointA by applying the SDRE control with α=10-6 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Figure 19
Control torques from the formulations based on pointA by applying the SDRE control with α=10-6 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)


### 5.2. Libration Control for the Equations of Motion Based on PointC

The same orbit and orbital rate are considered in this subsection, and satelliteD is also fixedly located at the middle point of the tether. Thus, there are still only two ordinary differential equations taken into account, and each equation is considered with a control torque, which will be determined by applying the SDRE control method.Based on the formulations of the SDRE control, the two equations of motion are factorized as(51)Aθθ¨+Bθθ˙+Cθθ=Du,where A(θ), B(θ), and C(θ) are 2×2 matrices and θ refers to θ=θ1θ2T. Besides, the components of the matrices and the vector are written as(52)A11θ=a1L12,A12θ=A21θ=a3L1L2c12,A22θ=a2L22,B11θ=2a1L1L˙1,B12θ=-2a3L1L˙1c12,B21θ=2a3L2L˙1c12,B22θ=-2a2L2L˙1,C11θ=μmARA-31-ρARCL1s1-1-ρAρBL1L2s1c2+μmBRB-3-ρARCL1s1-ρA1-ρBL1L2s1c2+μmDRD-3-ρARCL1s1+ρAρBL1L2s1c2+2a1L1L˙1ψ˙C+a3L1L¨1s1c2+L1L2ψ˙C+θ˙22s1c2θ1-1,C12θ=μmARA-31-ρAρBL1L2c1s2+μmBRB-3ρA1-ρBL1L2c1s2-μmDRD-3ρAρBL1L2c1s2+a3-L1L¨1c1s2-L1L2ψ˙C+θ˙22c1s2-2L1L˙1ψ˙Cc12θ2-1,C21θ=μmARA-3ρBRCL2s2+1-ρAρBL1L2s1c2+μmBRB-3-1-ρBRCL2s2+ρA1-ρBL1L2s1c2+μmDRD-3ρBRCL2s2-ρAρBL1L2s1c2+a3L2L¨1s1c2-L1L2ψ˙C+θ˙12s1c2+2L2L˙1ψ˙Cc12θ1-1,C22θ=-μmARA-31-ρAρBL1L2c1s2-μmBRB-3ρA1-ρBL1L2c1s2+μmDRD-3ρAρBL1L2c1s2-2a2L2L˙1ψ˙C+a3-L2L¨1c1s2+L1L2ψ˙C+θ˙12c1s2θ2-1.In order to compare the simulations with those in previous subsection, the initial conditions are the same, so the initial libration angles are given asθ1(0)=0.9649° and θ2(0)=-1.0351°. Based on the initial condition, the orbital radius is given as RC = 47,521 km, and the orbital rate is the same as ψ˙C=7.2923×10-5 rad/s. Besides, the parameters are given as R˙C=R¨C=ψC(0)=ψ¨C=0. In the SDRE control, the weighting matrices are the same as Q=diag⁡(1,1,0,0) and R=diag⁡(1,1). Note that the weighting matrices have different meaning due to different state variables. Thus, although most of parameters are the same as those in previous subsection, the simulation results will be different. Similarly, two cases are tested to demonstrate the prescribed degree of stability.Case 3 (α=0).
Figure20 shows the libration angles as functions of time, and the results show that the SDRE control can suppress the librations around within 10 orbits. Although the time responses seemingly have high-frequency oscillations in the transient states (see Figure 20(a)), the oscillation period is around 3,954.8 seconds (see Figure 20(b)), which is almost twice time in Figure 16(b). Figure 21 shows the control torques as function of time, and the results show that the maximum control torque is around 0.35 MN-m, which is also larger than the torque in Figure 16.Figure 20
Libration angles from the formulations based on pointC by applying the SDRE control with α=0 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Figure 21
Control torques from the formulations based on pointC by applying the SDRE control with α=0 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Case 4 (α=10-6).
Seemingly, the convergence rate forα=0 is little slow, so this test case is based on α=10-6. Figure 22 shows the libration angles as functions of time, and the results show that the SDRE control can suppress the librations around within 8 orbits, which indicates that the proposed SDRE control with a prescribed degree of stability can effectively reduce the converge time. This figure also shows a high-frequency oscillation in the transient state, and the oscillation period is the same as 3,954.8 seconds. Figure 23 shows the control torques as function of time, and the results show that the profile is similar to that in Case 1, and the difference is that the torque converges faster.Figure 22
Libration angles from the formulations based on pointC by applying the SDRE control with α=10-6 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Figure 23
Control torques from the formulations based on pointC by applying the SDRE control with α=10-6 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)


## 5.1. Libration Control for the Equations of Motion Based on PointA

One considers that the tethered satellite system moves on a circular orbit with a constant rotating rate and satelliteD is fixedly located at the middle point of the tether. Thus, there are only two ordinary differential equations taken into account, and each equation is considered with a control torque, which is obtained based on the SDRE control method.Based on the formulations of the SDRE control, the two equations of motion are factorized as(49)Aϕϕ¨+Bϕϕ˙+Cϕϕ=Du,where A(ϕ), B(ϕ), and C(ϕ) are 2×2 matrices, D is a 2×2 identity matrix, ϕ=ϕ1ϕ2T, and u=TATDT. Besides, the components of the matrices and the vector are written as(50)A11ϕ=mB+mDL12,A12ϕ=A21ϕ=mBL1L2c12,A22ϕ=mBL22,B11ϕ=2mB+mDL1L˙1,B12ϕ=-2mDL1L˙1c12,B21ϕ=2mDL2L˙1c12,B22ϕ=-2mBL2L˙1,C11ϕ=μmBRB-3-RAL1s1+L1L2s1c2+μmDRD-3-RAL1s1+mB+mD2L1L˙1ψ˙A+RAL1ψ˙A2-L1R¨As1+mBL1L2ψ˙A+ϕ˙22s1c2+L1L¨1s1c2ϕ1-1,C12ϕ=μmBRB-3-L1L2c1s2+mB+mD2L1R˙Aψ˙Ac1+L12+RAL1c1ψ¨A+mB-2L1L˙1ψ˙Ac12-L1L2ψ˙A+ϕ˙22c1s2+L1L2ψ¨Ac12-L1L¨1c1s2ϕ2-1,C21ϕ=μmBRB-3L1L2s1c2+mB-2L2L˙1ψ˙Ac12+2L2R˙Aψ˙Ac2-L1L2ψ˙A+ϕ˙12s1c2+L22+L1L2c12+RAL2c2ψ¨A+L2L¨1s1c2ϕ1-1,C22ϕ=μmBRB-3RAL1s2-L12c1s2+mB-2L2L˙1ψ˙A+RAL2ψ˙A2s2+L1L2ψ˙A+ϕ˙12c1s2-L2R¨As2-L2L¨1c1s2ϕ2-1.The simulations are performed based on the following parameters:RA = 42,522 km, ψ˙A=7.2923×10-5 rad/s, and R˙A=R¨A=ψA(0)=ψ¨A=0. Besides, the initial conditions are given as ϕ1=1°, ϕ2=-1°, and ϕ˙1=ϕ˙2=0. For simplicity, the weighting matrices in the SDRE control are selected as Q=diag⁡(1,1,0,0) and R=diag⁡(1,1). In order to demonstrate the performance of the prescribed degree of stability in the SDRE control, there are two cases tested as follows.Case 1 (α=0).
Figure16 shows the libration angles as functions of time, and the results show that the SDRE control can suppress the librations around within 10 orbits. Although the time responses seemingly have high-frequency oscillations in the transient states (see Figure 16(a)), the oscillation period is around 1,930 seconds (see Figure 16(b)). Figure 17 shows the control torques as function of time, and the results show that the maximum control torque is around 0.16 MN-m.Figure 16
Libration angles from the formulations based on pointA by applying the SDRE control with α=0 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Figure 17
Control torques from the formulations based on pointA by applying the SDRE control with α=0 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Case 2 (α=10-6).
Seemingly, the convergence rate forα=0 is little slow, so this test case is based on α=10-6. Figure 18 shows the libration angles as functions of time, and the results show that the SDRE control can suppress the librations around within 8 orbits, which indicates that the proposed SDRE control with a prescribed degree of stability can effectively reduce the converge time. This figure also shows a high-frequency oscillation in the transient state, and the oscillation period is the same as 1,930 seconds. Figure 19 shows the control torques as function of time, and the results show that the maximum control torque is around 0.28 MN-m, which is larger than those presented in previous case. This is the penalty to reduce the convergence time.Figure 18
Libration angles from the formulations based on pointA by applying the SDRE control with α=10-6 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Figure 19
Control torques from the formulations based on pointA by applying the SDRE control with α=10-6 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)


## 5.2. Libration Control for the Equations of Motion Based on PointC

The same orbit and orbital rate are considered in this subsection, and satelliteD is also fixedly located at the middle point of the tether. Thus, there are still only two ordinary differential equations taken into account, and each equation is considered with a control torque, which will be determined by applying the SDRE control method.Based on the formulations of the SDRE control, the two equations of motion are factorized as(51)Aθθ¨+Bθθ˙+Cθθ=Du,where A(θ), B(θ), and C(θ) are 2×2 matrices and θ refers to θ=θ1θ2T. Besides, the components of the matrices and the vector are written as(52)A11θ=a1L12,A12θ=A21θ=a3L1L2c12,A22θ=a2L22,B11θ=2a1L1L˙1,B12θ=-2a3L1L˙1c12,B21θ=2a3L2L˙1c12,B22θ=-2a2L2L˙1,C11θ=μmARA-31-ρARCL1s1-1-ρAρBL1L2s1c2+μmBRB-3-ρARCL1s1-ρA1-ρBL1L2s1c2+μmDRD-3-ρARCL1s1+ρAρBL1L2s1c2+2a1L1L˙1ψ˙C+a3L1L¨1s1c2+L1L2ψ˙C+θ˙22s1c2θ1-1,C12θ=μmARA-31-ρAρBL1L2c1s2+μmBRB-3ρA1-ρBL1L2c1s2-μmDRD-3ρAρBL1L2c1s2+a3-L1L¨1c1s2-L1L2ψ˙C+θ˙22c1s2-2L1L˙1ψ˙Cc12θ2-1,C21θ=μmARA-3ρBRCL2s2+1-ρAρBL1L2s1c2+μmBRB-3-1-ρBRCL2s2+ρA1-ρBL1L2s1c2+μmDRD-3ρBRCL2s2-ρAρBL1L2s1c2+a3L2L¨1s1c2-L1L2ψ˙C+θ˙12s1c2+2L2L˙1ψ˙Cc12θ1-1,C22θ=-μmARA-31-ρAρBL1L2c1s2-μmBRB-3ρA1-ρBL1L2c1s2+μmDRD-3ρAρBL1L2c1s2-2a2L2L˙1ψ˙C+a3-L2L¨1c1s2+L1L2ψ˙C+θ˙12c1s2θ2-1.In order to compare the simulations with those in previous subsection, the initial conditions are the same, so the initial libration angles are given asθ1(0)=0.9649° and θ2(0)=-1.0351°. Based on the initial condition, the orbital radius is given as RC = 47,521 km, and the orbital rate is the same as ψ˙C=7.2923×10-5 rad/s. Besides, the parameters are given as R˙C=R¨C=ψC(0)=ψ¨C=0. In the SDRE control, the weighting matrices are the same as Q=diag⁡(1,1,0,0) and R=diag⁡(1,1). Note that the weighting matrices have different meaning due to different state variables. Thus, although most of parameters are the same as those in previous subsection, the simulation results will be different. Similarly, two cases are tested to demonstrate the prescribed degree of stability.Case 3 (α=0).
Figure20 shows the libration angles as functions of time, and the results show that the SDRE control can suppress the librations around within 10 orbits. Although the time responses seemingly have high-frequency oscillations in the transient states (see Figure 20(a)), the oscillation period is around 3,954.8 seconds (see Figure 20(b)), which is almost twice time in Figure 16(b). Figure 21 shows the control torques as function of time, and the results show that the maximum control torque is around 0.35 MN-m, which is also larger than the torque in Figure 16.Figure 20
Libration angles from the formulations based on pointC by applying the SDRE control with α=0 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Figure 21
Control torques from the formulations based on pointC by applying the SDRE control with α=0 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Case 4 (α=10-6).
Seemingly, the convergence rate forα=0 is little slow, so this test case is based on α=10-6. Figure 22 shows the libration angles as functions of time, and the results show that the SDRE control can suppress the librations around within 8 orbits, which indicates that the proposed SDRE control with a prescribed degree of stability can effectively reduce the converge time. This figure also shows a high-frequency oscillation in the transient state, and the oscillation period is the same as 3,954.8 seconds. Figure 23 shows the control torques as function of time, and the results show that the profile is similar to that in Case 1, and the difference is that the torque converges faster.Figure 22
Libration angles from the formulations based on pointC by applying the SDRE control with α=10-6 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)Figure 23
Control torques from the formulations based on pointC by applying the SDRE control with α=10-6 for (a) the first 10 orbits and (b) the first 0.1 orbits.
(a) (b)


## 6. Conclusions

This paper investigates the nonlinear dynamic analysis and control of a tethered satellite system. First, the system is modeled by a two-piece dumbbell. There are three satellites in the system, and each satellite is assumed to be a point mass. A constant-length massless tether connects the three satellites. Based on the model, two sets of equations of motion are derived, which are based on the first satellite and the mass center of the system. Secondly, the SDRE with a prescribed degree of stability is proposed, which can reduce the convergence time of a nonlinear system. Thirdly, the nonlinear dynamic analysis is performed by numerical simulations. The results show that the librations of the system affected the motion of the middle satellite. Also, its motion affects the orbital radius and orbital rate. Finally, the proposed SDRE approach is developed and is applied to the system, and the numerical examples demonstrate the efficiency of the proposed approach in the system.



---

*Source: 3610507-2016-09-04.xml*