Nonlinear dynamics and reconfiguration control of two-satellite Coulomb tether formation at libration points

JID:AESCTE AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.1 (1-12)

Aerospace Science and Technology ••• (••••) •••–•••

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Nonlinear dynamics and reconfiguration control of two-satellite Coulomb tether formation at libration points Jing Huang ∗ , Guangfu Ma, Gang Liu School of Astronautics, Harbin Institute of Technology, Harbin 150001, China

a r t i c l e

i n f o

Article history: Received 8 February 2014 Received in revised form 10 June 2014 Accepted 18 June 2014 Available online xxxx Keywords: Coulomb tether formation Libration points Restricted three-body problem Nonlinear robust control

a b s t r a c t The use of inter-satellite electrostatic (Coulomb) forces for formation flying is attractive due to the prominent advantages of no propellant consumption or plume impingement issues. This paper investigates the nonlinear attitude–orbit coupling dynamics and control for the reconfiguration of a twosatellite Coulomb tether formation near Earth–Moon libration points. The nonlinear dynamic model is derived by utilizing analytical mechanics theory, and analysis on environment disturbance related to the solar radiation pressure effects is put forward. Considering the high nonlinearity and uncertainty in the dynamic model, a closed-loop control strategy based on the indirect robust control scheme and the θ –D method is proposed with precise and good robust performance. Here, the Coulomb force is used for the longitudinal direction control, and the inertial micro-thrusters are activated for control in the transverse directions. Theoretical analysis and numerical simulation results are presented to validate the feasibility of the proposed dynamics model and proposed control strategy for the robust reconfiguration and station-keeping mission. © 2014 Elsevier Masson SAS. All rights reserved.

1. Introduction In the last few years, the concept of using electrostatic Coulomb forces to control satellite formation flight of close-proximity is a promising form of propulsion that exploits the electric potential to generate desired inter-satellite forces [15,14]. Certain benefits ranging from fuel consumption to high precision can result from a Coulomb tether formation. Cover [5] studied the active charge control as early as 1966, which employed the electrostatic forces to inflate and maintain the shape of a large reflecting structure. For tight formation control, such nearly-propellant-less thrust is advantageous over conventional electric propulsion regarding the thruster plume contamination of the neighboring satellite. Moreover, the Coulomb propulsion has some advantages such as fast throttling (transition time ∼ms) and remarkably little electrical power required, which enable high precision satellite formation flying with long lifetime. Some researches about the application of the Coulomb tether satellite system have been done in high accuracy optical interferometry missions, advanced docking/rendezvous, autonomous inspection, and the deployment/removal of

*

Corresponding author at: School of Astronautics, Harbin Institute of Technology, Harbin, 150001, China. Currently, a visiting scholar at Department of Mechanical and Industrial Engineering, Concordia University, Montreal, Quebec H3G 1M8, Canada. E-mail address: [email protected] (J. Huang). http://dx.doi.org/10.1016/j.ast.2014.06.001 1270-9638/© 2014 Elsevier Masson SAS. All rights reserved.

instrument [2], etc. However, there are some inherent drawbacks of this propulsion system. The formation dynamics are nonlinear and highly coupled, and the internal Coulomb forces depend on every satellite’s position and charge. For this reason, Coulomb satellite system requires a careful balance between the inter-satellite forces and the orbital dynamics. Also, as the electrostatic forces are internal to the formation, external forces such as thrusters are needed to reorient the whole Coulomb formation to a new orientation. There have already been a number of researches about the dynamical modeling and control of conventional tethered satellite system [4]. For the control of Coulomb tether formation, Parker and King [15,14] first introduced the novel method of exploiting Coulomb forces to control free-flying spacecrafts in 2002. The analytic open-loop solutions were obtained for Hill-frame invariant static Coulomb formations with symmetry assumptions. All relative motions of the charged spacecrafts were eliminated because of the Coulomb forces. Based on their pioneering work, there have been some investigations on the dynamics and control of the Coulomb system. In Refs. [2,27], the systematic analytic solutions for some Coulomb satellite formations were presented. Besides, Berryman [1] numerically demonstrated the possibility of the spacecraft formation with as many as nine satellites using Coulomb forces in GEO. These open-loop solutions for static Coulomb formations are all dynamically unstable without a feedback control law to stabilize the motion. To solve the problem, Natarajan and Schaub [23]

JID:AESCTE

AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.2 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

2

Nomenclature ao Ai C Cr F Fs F sr F srd F sd Gp g max (·) Is I¯ C i, j, k kc Ld L d1 , L d2 L l1 , l2 M1 M2 m

acceleration of the origin O cross section area of the satellite i that is facing the Sun center of mass of system radiation pressure coefficient solar radiation flux inertial force vector related to solar radiation pressure solar radiation pressure force differential solar radiation pressure effect force vector differential inertial force vector related to the solar radiation pressure common centroid in a circular orbit upper limit of uncertain disturbance moment of inertia about a line normal to the virtual tether by C of Coulomb satellite system inertial tensor of the system at C the unit vectors along the X -, Y -, and Z -axes Coulomb’s constant distance between the two satellites distances of each satellite from C distance between the primaries the distances of M 1 and M 2 from G p mass of the first primary (Earth) mass of the second primary (Moon) total mass of the Coulomb satellites system

presented closed-loop feedback control laws to stabilize the virtual Coulomb tether system for in-orbit two-craft configurations. Based on the linearized time-varying dynamical models, they investigated the linear dynamics and stability of the system with the proposed controller for all three equilibrium configurations. However, the linearized model is not accurate enough in performing reconfigurations in the nonlinear regime. For the control of the nonlinear model, Natarajan [22] presented a feedback control law to reconfigure the two-craft Hill-frame Coulomb formation. Then, Inampudi [10] applied the pseudo-spectral discretization algorithm and the numerical optimization method to minimize the time and fuel consumption for two-craft Coulomb formation reconfigurations in the circular orbits. Huang [9] developed the generalized dynamics of multi-spacecraft formation in a central gravitational field with non-contacting internal forces, such as Coulomb forces, based on Kane method. The necessary condition and constraints to achieve static formation were derived by examining the relative equilibrium. Jones [13] introduced invariant manifold theory to two-craft Coulomb formation configurations in the two-body gravity Hill-frame model. In his research, the original control problem was converted to a parameter optimization problem, which was numerically solved through particle swarm optimization. The dynamics and control of a tether satellite system in a three-body gravitational environment have received a lesser attention than the two-body case. A location near the libration points is especially suitable for remote sensing mission. It is feasible to exploit the stability condition near L 2 point to further reduce the fuel consumption. Pioneering studies of conventional tethered satellite system in the deep space were done by Gates [7] and later by [6]. These studies investigated the dynamics model of the multi-tethered systems in arbitrary configurations near collinear libration points without considering the gravitational attractions. Subsequent works were done by Misra [20] and Wong [28], but the analysis they did was only on the planar motion of the system near equilibrium positions. Recently, Sanjurjo-Rivo, Lucas, and Peláez [26] investigated the dynamics of a tethered satellite sys-

mi mass of the satellite i, i = 1, 2 P n (·) Legendre polynomials qi satellite charge Q˜ positive semi-definite weighting matrix ˜ R positive definite weighting matrix R 1 , R 2 position vectors of C with respect to M 1 and M 2 R dm1 , R dm2 the scalar distances of dm from M 1 and M 2 rs position vector from the Sun to the satellite in AU s direction of the incident sun vector Tc control force vector Td disturbance vector TL Coulomb force u 1 , u 2 , u 3 unit vectors of the body frame attached to the virtual tether vC velocity vector of C X L2 X position of L 2 point in the frame OXYZ α, β in-plane and out-plane libration angles γ,δ cone and clock angles of the satellite-normal relative to the orbit frame λd Debye length ω angular velocity of revolution of the two primaries ωB angular velocity of Earth–Moon system around the Sun ωoxyz angular velocity vector of the frame OXYZ Ω angular velocity vector of the Coulomb satellite system relative to the frame OXYZ

tem near collinear libration points. A simple feedback control law was proposed to achieve system stabilization. Cai [3] developed the dynamical formulation of a rotating triangular tethered satellite formation near libration points and performed the parametric studies of some effects to the system during the deployment and retrieval stages. However, the Hill approximation approach they used is a simplification of the Circular Restricted Three Body Problem (CRTBP). For better analysis on the real system dynamics and control, a high-precision description as the CRTBP is needed. Liu [18] investigated the dynamics of a rotating two-body tethered satellite system near the collinear libration points in the CRTBP. A θ –D nonlinear suboptimal output tracking controller is developed for the stabilization mission of the system along the Halo orbits. However, they did not consider any disturbances on the system. For the Coulomb virtual tether satellite system in the deep space, the NIAC report [14] first analyzed the feasibility of Coulomb force control for a static collinear five-vehicle formation in the vicinity of Sun–Earth libration points. Pettazzi [25] investigated a distribute navigation technique called Equilibrium Shaping, which was designed to drive a swarm of Coulomb satellites with electric propulsion to a desired configuration at the libration points. However, the control algorithm does not have analytical stability guarantees. The Debye length, which relates to the lower bound on the electrostatic field strength of plasma, varies between 10–40 m [8] when the spacecraft is moving around the collinear Earth–Moon libration points. It constrains the maximum possible formation length. Although the formation length is relatively small compared with the conventional tether, multi-spacecraft equilibrium formations still exist at these points. Besides, charged relative motion dynamics, which are primarily influenced through classical electrostatic, cannot be ignored [21]. These results motivate us to study the dynamics and control of a two-craft Coulomb formation at Earth–Moon libration points. To solve this problem, Inampudi [11] concluded that three equilibrium configurations exist for a two-craft Coulomb formation at

JID:AESCTE AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.3 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

the gravitational three-body libration points. Then for such formation at Earth–Moon collinear libration points, a generic nonlinear controller was proposed to counteract the effects of the solar perturbation in Ref. [12]. It should be noted that in Refs. [11,12], the center of mass of the Coulomb tether satellite system considered was fixed at the collinear libration points, while the interaction effects of the attitude, distance between the satellites and orbit motion were not considered. It is necessary to consider the motion of the system position due to the effects of coupling dynamics, gravity and disturbance. The focus of the paper is to maneuver the two-satellite formation between two configurations and keep the system staying steadily at Earth–Moon L2 point with considering the system orbit motion, which is different from Inampudi’s work, and also the major contribution of the paper. The motion of a Coulomb tether satellite system can be described by the equations of the conventional tethered satellite system in the CRTBP, with the addition of the inter-craft force produced by the Coulomb propulsion. Differential solar drag is the largest disturbance acting on a tether formation at libration points. To deal with the external disturbance, a transformation method inspired by the indirect robust control presented in Refs. [17,30] is proposed. The uncertain dynamic system is changed to a nominal system of optimal control problem with a modified cost function including uncertainty bounds. However, in most cases, it is impossible to obtain an analytical solution for general nonlinear optimal control problem. Therefore, a control approach based on the θ –D method [29] is chosen for the transformed nonlinear optimal control problem. It is improved to solve an optimal tracking control problem with uncertainty. Moreover, the differential solar drag disturbance is analyzed and transferred to a quadratic form such that the robust optimal control method can be applied. The remainder of this paper is organized as follows. In Section 2, the attitude and orbit coupling dynamics of a two-satellite Coulomb tether system in the CRTBP is investigated. After that, the environment disturbance related to the solar radiation pressure effects is discussed. In Section 3, an indirect robust control design scheme is introduced, and both robust stability and optimality performance can be achieved. By combining the indirect robust control scheme with the θ –D method, a closed-form nonlinear robust suboptimal controller is obtained to solve the Coulomb satellite formation reconfiguration control and station-keeping problem. At last, the performance of this control scheme is evaluated through numerical simulation in Section 4 and the conclusion is given in Section 5. 2. System model 2.1. Description of the system In this paper, a new approach for the general treatment of a two-satellite Coulomb tether system near the Earth–Moon collinear libration points is considered. To proceed, throughout this work, the following assumptions are made. (1) The inter-satellite force undergoes both tensile and compressive forces along the line-of-sight direction between the two satellites. (2) The satellites are considered to be point masses moving around the libration points. The gravitational attraction between the two satellite masses is neglected. (3) The dynamic model is based on the CRTBP with the Earth as one primary and the moon as the second, the two primary revolve around their common centroid G p in a circular orbit.

3

Fig. 1. Geometry of the Coulomb tether system in the CRTBP.

The Coulomb satellite system under consideration consists of two end masses, denoted by m1 and m2 , and the distance between each other is L d . The center of mass of the system is defined as C . As shown in Fig. 1, the two-satellite Coulomb system is seen as the smallest body in the CRTBP, the two primary masses are M 1 and M 2 . l denotes the distance between the two primaries, while l1 and l2 denote the distances of M 1 and M 2 from G p , respectively. It can be calculated that l1 = ν l and l1 = (1 − ν )l, where ν = M 2 /( M 1 + M 2 ). The motions of the Coulomb satellite system are described using two sets of rotating coordinate axes. The primary frame [G p XYZ] has its origin at G p and X -axis points from M 1 to M 2 . Y -axis is perpendicular to X -axis in the plane of motion of the primary bodies. The synodic frame [G p XYZ] is rotating around the axis G p Z with a constant angular rate ω , with ω =



G ( M 1 + M 2 )/l3 , G is the universal gravitational constant. Since the motion of Coulomb satellite system in the vicinity of M 2 is mainly studied in this paper, a secondary set of axes, [OXYZ], is located at the center of mass of the secondary primary of interest and is parallel to the corresponding axis of the frame [G p XYZ]. The unit vectors along the X -, Y -, and Z -axes are denoted by i, j, and k, respectively. The position vectors of C relative to M 1 and O , respectively, can be written as

R 1 = ( X + l)i + Y j + Z k R2 = X i + Y j + Z k

R1 = | R 1| =

R2 = | R 2| =





( X + l)2 + Y 2 + Z 2 (1)

X2

+

Y2

+

Z2

(2)

2.2. Equations of motion The equations of motion of the two-satellite Coulomb tether system with hybrid thrusting (both electrostatic and inertial thrusting) are derived based on the method of Lagrangian. For the motion relative to the synodic frame OXYZ, the kinetic energy of the system is

K=

=

1 2 1 2

mv 2C + mv 2C

+

1 2 1 2

1

1

2 2 Ω ◦ I¯ C ◦ Ω + m1 L˙ d1 + m2 L˙ d2

2

I s u˙ 21

+

1 2

2 m1 L˙ d1

2

+

1 2

2 m2 L˙ d2

(3)

In this expression, Ω is the angular velocity of the Coulomb satellite system relative to the frame OXYZ; I¯ C is the inertia tensor of the system at C ; v C is the velocity vector of C and v C is the module of v C ; L d1 and L d2 are the distances of each satellite to the center of mass C , respectively. I s is the moment of inertia about a line normal to the virtual line connecting the two satellites by the center of mass C , given by I s = mL d2 b, where b = sin2 2φ/4, φ

is the mass angle defined by sin2 φ = m2 /m, cos2 φ = m1 /m [24]. As shown in Fig. 2, the attitude dynamics of the system is described using the libration angles (α , β ) as generalized coordinates.

JID:AESCTE

AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.4 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

4

V g2 = −

∞ G M1 

R1

(−1)n P n (cos χ1 ) I n1

n =0

∞ G M2 

−

R2



(−1)n P n (cos χ2 ) I n2

(11)

n =0

ηni dm, i = 1, 2

I ni =

(12)

m Fig. 2. Frame C u 1 u 2 u 3 attached to the tether when attitude is described with libration angles α (in-plane) and β (out-plane): (a) Definition of tether angles and u 1 . (b) Definition of vector u 2 on X Y plane. (c) Definition of vector u 3 on N Z plane.

In teams of the libration angles, the body frame defined by the three unit vectors (u 1 , u 2 , u 3 ) can be expressed as follows

u 1 = cos β cos α i + cos β sin α j + sin β k

I 0i = m,

V g2 ≈ −

 1   1  2 ˙ 2 cos2 β m X˙ + Y˙ 2 + Z˙ 2 + mL d2 b β˙ 2 + α 2 2 1 + m L˙ d2 sin2 2φ 8

V c = kc (5)

where ao is the acceleration of the origin O , ωoxyz is the angular velocity of the frame OXYZ. Then the following expression of generalized potential can be derived



2

2

V g1 = − mω X + (1 − ν )l + Y 2 1 − I s ω2 cos2 β − I s ωα˙ cos2 β 2



V g2 = − m

G M1 R dm1



dm − m

G M2 R dm2

dm

− mω( X Y˙ − Y X˙ ) (7)

(8)

where R dm1 and R dm2 are the scalar distances of dm from the attracting centers M 1 and M 2 , respectively. By introducing the variables η1 = s1 / R 1 , η2 = s2 / R 2 , s1 and s2 are the distances of the mass element dm from the mass center of the body, in a general case, ηi is small: ηi  1, i = 1, 2. As a consequence,

1 R dm1

=

R 1 1 + 2η1 cos χ1 + η12

= 1 R dm2

=

1



∞ 1 

R1

(−1)n 1n P n (cos

η

χ1 )

(9)

n =0

R 2 1 + 2η2 cos χ1 + η22

=

∞ 1 

R2

(−1)n 2n P n (cos

η

G M 2m



2

R1

1+

R2

d ∂ L(q, q˙ )

∂ q˙

Ld

 b P 2 (cos χ1 )

2

 b P 2 (cos χ2 )

R2

(14)

(10)

n =0

where χi , i = 1, 2 is the angle between the vectors R i and u 1 , cos χ1 = ( R 1 · u 1 )/ R 1 , cos χ2 = ( R 2 · u 1 )/ R 2 , P n (·) are the Legendre polynomials. Eq. (8) can be expressed as the following series

(15)

−

∂ L(q, q˙ ) =Q ∂q

(16)

where q = [ X , Y , Z , α , β, L d ] T is a generalized coordinate, Q = [ Q X , Q Y , Q Z , Q α , Q β , Q L ] T = T c + T d is the corresponding generalized force vector including the external control force vector T c = [ T X , T Y , T Z , T α , T β , 0] T and disturbance vector T d = [d X , d Y , d Z , dα , dβ , d L ] T . To prevent numerical difficulties with small quantities, a set of dimensionless quantities are defined as follows

x=

X l

y=

,

Y l

,

z=

Z l

ε=

,

Ld l

,

τ = ωt

(17)

From there on, the time derivative q stands for the derivative with respect to the non-dimensional time τ : q = dq/dτ . Making use of Eqs. (16) and (17), the equations of motion of two-satellite Coulomb tether formation for the CRTBP in the dimensionless form can be expressed as

x − 2 y  − (1 − ν + x) +



+

χ2 )

1+

Ld

q1 q2 − λLd e d Ld

= εb

1



where q i , i = 1, 2 is the satellite charge, and the parameter kc = 8.99 × 109 Nm2 /C2 is Coulomb’s constant. The exponential term in the Coulomb potential depends on the Debye length parameter λd which controls the electrostatic field strength of plasma shielding between the satellites. The nonlinear equations of motion are deduced from the Lagrangian L = K − ( V g1 + V g2 + V c ) of the system in the following form

dt

The potential of gravitational forces is classically given by



(13)

where P 2 (x) is given by P 2 (x) = (3x2 − 1)/2. The Coulomb potential for the two-satellite formation is

  1 V g1 = − mR 2 2ao + ωoxyz × (ωoxyz × R 2 ) + 2ωoxyz × v C 2 1 − ωoxyz ◦ I¯ C ◦ (ωoxyz + 2Ω) (6) 2

2



R1

(4)

The generalized potential takes the form

1

G M 1m

−

This way we obtain the following expression

K=

i = 1, 2

After the approximations for I ni , the potential of gravitational forces finally becomes

u 2 = − sin α i + cos α j u 3 = − sin β cos α i − sin β sin α j + cos β k

I 2i = I s / R 2i ,

I 1i = 0,

r13

+

xν r23

(1 − ν )  r15



3(r 1 · u )(i · u 1 ) − (1 + x) S 2 (cos χ1 )

  QX 3(r 2 · u 1 )(i · u 1 ) − xS 2 (cos χ2 ) + 2

ν r25

mω l

y  + 2x − y +

 = εb

(1 + x)(1 − ν )

y (1 − ν ) r13

(1 − ν )  r15

+

yν r23



3(r 1 · u 1 )( j · u 1 ) − y S 2 (cos χ1 )

(18)

JID:AESCTE AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.5 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••



ν

+



3(r 2 · u 1 )( j · u 1 ) − y S 2 (cos χ2 )

r25

z(1 − ν )

z +

r13

 = εb

(19)

r23

(1 − ν ) 



3(r 1 · u 1 )(k · u 1 ) − zS 2 (cos χ1 )

r15



3(r 2 · u 1 )(k · u 1 ) − zS 2 (cos χ2 )

r25

QY mω 2 l

zν

+

ν

+

+



 +



QZ mω 2 l

(20)

α  cos β − 2 1 + α  β  sin β =−

I s  Is

+3

 

1 + α cos β + 3

ν r25

(1 − ν ) r15

(r 2 · u 1 )(r 2 · u 2 ) +

(r 1 · u 1 )(r 1 · u 2 ) Qα

(21)

ε2 bmω2l2 cos β

Qβ

+ 2 ε bmω2l2  2  2   ε − ε β + 1 + α  cos2 β 2ε (1 − ν ) P 2 (cos χ1 ) 2εν P 2 (cos χ2 ) − − − 3 3 r1

r2

QL

(22)

TL bmω2l



S 2 (x) =

2

2

5x − 1



(24)

Here, the Coulomb forces are used for the longitudinal direction control, and the inertial micro-thrusters are activated for control in the transverse directions. The variable T L associated with Coulomb propulsion is expressed as

TL =

kc q1 q2 − λLd e d L d2

1+

Ld



(25)

λd

Since an infinite number of solutions satisfying Eq. (25) will be found, charges of the same magnitude across the satellite are chosen in this paper. For instance, for an equilibrium configuration assuming equal charges in magnitude and using Eq. (25) yields

q1 q2 = q1 =



T L L d2 λd

Ld

kc (λd + L d )

e λd



|q1 q2 |,

(26)

q2 = q1 ,

if T L ≥ 0

q2 = −q1 ,

if T L < 0

T

(28)

The inertial force vector F s in N due to the effects of solar radiation pressure (SRP) is expressed as [19]

where r 1 = (x + 1)i + y j +zk, r1 = |r 1 | = (x + 1)2 + y 2 + z2 , r 2 = xi + y j + zk, r2 = |r 2 | = x2 + y 2 + z2 and the function S 2 (x) is given by

3

axis G p Z with a constant angular velocity ω . Simultaneously, the Earth–Moon system orbits the Sun with angular velocity ω B . As a result, the incident sun line rotates in the orbit frame with the angular velocity ωs = ω − ω B . The direction of the incident sun vector s will vary continuously with respect to the synodic frame [G p XYZ] as

s = cos(ωs t ), − sin(ωs t ), 0

(23)

bmω2l

Fig. 3. Cone and clock angles (δ, γ ) of the satellite-normal relative to the orbit frame.



2 I s   β + 1 + α  sin β cos β Is (1 − ν ) ν (r 1 · u 1 )(r 1 · u 3 ) + 3 5 (r 2 · u 1 )(r 2 · u 3 ) =3 5 r1 r2

β  +

=

5

(27)

2.3. Environment disturbances In this section, the external disturbances of the Coulomb satellite system are discussed. Differential solar drag is the largest disturbance acting on a tether formation at libration points [19]. Compared with the solar radiation pressure effects, the other environmental disturbances become relatively insignificant and can be neglected. In the vicinity of the collinear libration points of the Earth–Moon system, the sun lines are treated as parallel lines. As discussed earlier, the synodic frame [G p XYZ] is rotating around the

Fs = −

Cr A F

rs

c

r s 3

(29)

where r s is the position vector from the Sun to the satellite in AU, A is the cross section area of the satellite that is facing the Sun in m2 . The constant F = 1372.5398 Watts/m2 is the solar radiation flux, c = 2.997 × 108 m/s is the speed of light, and C r = 1.3 is the radiation pressure coefficient. The components of the SRP force F sr = [ F sr1 , F sr2 , F sr3 ] T in the Earth–Moon rotating frame are given by

F sr1 = F s cos2 γ cos(ωs t − γ ) F sr2 = − F s cos2 γ sin(ωs t − γ ) sin δ F sr3 = F s cos2 γ sin(ωs t − γ ) cos δ

(30)

where F s = F s . F sr1 , F sr2 and F sr3 are the components in orbit radial direction, orbital velocity direction and the orbit normal direction, respectively. The cone angle and the clock angle, which are denoted by δ and γ , indicate the orientation of each satellite with respect to the orbit frame (see Fig. 3). In this study, these angles for each satellite are fixed. It can be seen that the SRP force in the Earth–Moon system is periodic and time varying. Let us introduce the differential solar radiation pressure effect force vector F srd



F srd = F sd cos2 γ cos(ωs t − γ ),

− F sd cos2 γ sin(ωs t − γ ) sin δ, T F sd cos2 γ sin(ωs t − γ ) cos δ

(31)

where F sd is the differential inertial force vector related to the solar radiation pressure and F sd = F sd . Then, the non-dimensional disturbance d due to the differential solar radiation pressure effects can be derived from Eqs. (30) and (31)



d = 01×6 ,

F sr1

,

F sr2

,

F sr3

mω 2 l mω 2 l mω 2 l

F srd · u 3 bmω2l

T

,

F srd · u 1

,

F srd · u 2

ε2 bmω2l2 cos β ε2 bmω2l2

, (32)

The disturbance will aggravate the instability of the Coulomb formation at the collinear libration points as well as affect the configuration of the formation. Therefore, a feasible control strategy is required to maintain the Coulomb system in the vicinity of libration points and maneuver the two-satellite formation between two configurations.

JID:AESCTE

AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.6 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

6

3. Suboptimal robust controller design for the reconfiguration of Coulomb formation In this section, an indirect method is designed to convert the nonlinear robust control problem to an equivalent optimal control problem. Then, the θ –D suboptimal control approach is developed to obtain a closed-form feedback control law for the reconfiguration mission of the Coulomb satellite system near libration points. 3.1. Problem statement The dynamic equations (18)–(23) can be rewritten in the following state-space form

X  = f ( X ) + B ( X )U ( X ) + d( X )

(33)

where the system variable X is defined as X = [x, y , z, α , β, ε , x ,



y  , z , α  , β  , ε  ] T , control input U = TL bmω2 l

T

, the external disturbance vector



d = 01×6 ,

dX

2 bm

dY

,

dZ

,

mω 2 l mω 2 l mω 2 l

dβ

ε

T TX , TY , T Z , Tα , β , mω2 l mω2 l mω2 l bmω2 bmω2

ω

2l2

,

ε2 bmω2l2 cos β

T

dL bmω

dα

,

⎢1 ⎢ 1 ⎢ ⎢ 1 ⎢ B(X ) = ⎢ ⎢ ⎢ ⎣

1 ε 2 l2 cos2 β

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

1

ε 2 l2

(34)

Proof. Consider the optimal control problem (36) and (37). Let U ∗ , v ∗ be the optimal solution, and define

∞ 

2 ρ 2 gmax ( X ) + e T Q˜ e

U ,v

τ

2

∗T

T

˜ R˜ + R





0

subject to

(43)

+

(44)

 I − B(X )B (X ) = 0

The optimal controller can be obtained from these three equations. The capability of the controller to stabilize the system (35) will be proved using Lyapunov’s asymptotic stability theorem. V ( X ) is continuously differentiable with respect to X . V ( X ) > 0 when X = 0 and V ( X ) = 0 when X = 0. Therefore, V ( X ) is chosen as a Lyapunov function for the original robust control problem. Evaluating the time derivation of V ( X ) and using the Eqs. (41)–(43), we obtain



    = V XT f ( X ) + B ( X )U ∗ + I − B ( X ) B + ( X ) v ∗   + d − V XT I − B ( X ) B + ( X ) v ∗

1 1 1  2 T 2 = − ρ 2 gmax ( X ) − e T Q˜ e − U ∗T R˜ R˜ U ∗ + ρ 2  v ∗  2

2

2 ∗T

−ρ v

d − 2U

∗T

2

˜ T R˜ B + ( X )d R

T

(45) T

˜ R˜ U ∗ − 2U ∗ T R˜ R˜ B + ( X )d, we have For the term −U ∗ T R T T −U ∗T R˜ R˜ U ∗ − 2U ∗T R˜ R˜ B + ( X )d  2  2 = 2 R˜ B + ( X )d −  R˜ U ∗ + B + ( X )d  2 ≤ 2 R˜ B + ( X )d

Thus, the time derivation of V ( X ) becomes

(37)

The error vector e ∈ R is n

(38)

(42)

=0



2 ρ 2 gmax ( X ) + e T Q˜ e + 2U T R˜ R˜ U + ρ 2 v 2 dτ (36)

 X = f ( X ) + B ( X )U + I − B ( X ) B ( X ) v

+

V XT

V XT B ( X )

V  = V XT X  = V XT f ( X ) + B ( X )U ∗ + d (35)

 

    + V XT f ( X ) + B ( X )U ∗ + I − B ( X ) B + ( X ) v ∗ = 0

ρ v

n

+



2 ∗T

T

(41)

T 2 2 ρ 2 gmax ( X ) + e T Q˜ e + 2U ∗T R˜ R˜ U ∗ + ρ 2  v ∗ 

2U

1

∞ 1 

e = X − Xr

where ξ12 ≤ ξ 2 = λmin ( Q˜ ), and λmin (·) means the minimum eigenvalue. According to Ref. [17], for the matched uncertainty, when the weighting matrices in performance index are constant scalars, the solution to optimal control problem is also a solution to robust control problem. This conclusion is established as above. Now the proof is given below.

1

where X ∈ R , U ( X ) ∈ R .d( X ) ∈ R is the uncertain disturbance with d( X ) ≤ g max ( X ). Define X r ∈ Rn as the desired state vector. It is to find the feedback control law U ( X ) that makes the state X in system (35) globally asymptotically converge to desired state X r . The above robust control problem can be solved indirectly by transforming it into an optimal control problem: minimize



(40)

as the optimal value function. From the optimal control theory [16], the following equations are obtained:

⎤

X  = f ( X ) + B ( X )U ( X ) + d( X )



2

 T + 2U R˜ R˜ U + ρ 2 v 2 dτ

The problem is defined as follows: Consider the nonlinear system described in the form of

2



T

06×6

m

(39)

2 ρ 2 v 2 + 2 R˜ B + ( X ) gmax ( X ) < ξ12 e 2

.

2l

BT

and ρ a positive tunable design parameter. The following condition is satisfied:

2

3.2. Indirect robust control design scheme

J=

− 1

 1 V X (τ ) = min

The uncertainty is assumed to be bounded as follows: d( X ) ≤ g max ( X ), where g max ( X ) is the upper limit of uncertain disturbance. Therefore, the problem addressed here is to find a feedback control law U ( X ) such that the states in system (33) will asymptotically converge to desired states with the bounded uncertainty d( X ).

n



B+ = B T B



,

The control input matrix B ( X ) is given by

⎡

˜ ∈ Rm×m positive where v is the auxiliary control, Q˜ ∈ Rn×n , R + definite matrices and B the pseudo-inverse of B given by

V ≤−

1

1

2 ρ 2 gmax ( X ) − e T Q˜ e

2 1

+ ρ 2





2  ∗ 2

v

2

 2 − ρ 2 v ∗T d + 2 R˜ B + ( X )d

(46)

JID:AESCTE AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.7 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 1 1  2 2 ≤ − ρ 2 gmax ( X ) − e T Q˜ e + ρ 2  v ∗  2 1

2

According to the condition (40), one can reduce the steady state error to achieve the requirement for the indirect control by selecting ˜ Then, the matrices appropriate weighting parameters ρ , Q˜ and R. Q and R in the cost function (54) are adjusted through the simulations to get a faster response while the control input is feasible. From the optimal control theory [16], the Hamiltonian of the above optimal control problem is

2

 2  2  + ρ 2  v ∗  + d 2 + 2 R˜ B + ( X )d 2 1

 2 1 2 ≤ − ρ 2 gmax ( X ) − e T Q˜ e + ρ 2  v ∗  2 1

2

 2 2 + ρ 2 gmax ( X ) + 2 R˜ B + ( X )d

2  2  2 1 T = − e Q˜ e + ρ 2  v ∗  + 2 R˜ B + ( X )d 2

1

H=

 2  2 ≤ −λmin ( Q˜ ) e 2 + ρ 2  v ∗  + 2 R˜ B + ( X )d

(47)

Since the condition (40) is satisfied, the following inequality can be derived





7

V  ≤ − ξ 2 − ξ12 e 2 < 0

2

eT Q e +

3.3. Suboptimal robust controller design for the Coulomb satellite system

2

M T R M + λT f ( X ) + λT B m ( X ) M

(55)

where λ ∈ R12 is the co-state vector and the optimal controller computed using the equation T M ∗ = − R −1 B m ( X )λ

(56)

The co-state equation takes the form of

(48)

According to the corollary of Barbalat Lemma, we can obtain limt →∞ e = 0. Then, the state in system (35) will achieve desired state globally. Therefore, the optimal control law U ∗ ( X ) obtained from the optimal control problem is also a robust control law which makes the closed-loop system globally asymptotically stable.

1

λ = −

T

∂H ∂ f (X) = −Q e − λ ∂X ∂X

(57)

The above nonlinear equation is extremely difficult to solve by direct methods. Hence, the θ − D technique is applied here, which provides an approximate closed-form solution by introducing instrumental variables to the cost function

J=

1

∞



e

2

T

Q +

∞ 



Diθ

i



T

e + M R M dτ

(58)

i =0

0

∞

The disturbance d( X ) is difficult to determine due to the change of the cross section area of the satellites facing the Sun. In this paper, the bound for d( X ) is chosen to be a quadratic function of the states such that the robust optimal control method can be applied

i where terms of i =0 D i θ , i = 1, 2 . . . n is a perturbation series in ∞ an auxiliary variable θ . D i and θ are chosen so that Q + i =0 D i θ i is positive semi-definite. Rewrite the state equation (53) in a linear factorization structure as follows [29]

  d( X )2 ≤ [x, y , z, α , β, ε ]σ [x, y , z, α , β, ε ]T

X  = f ( X ) + Bm ( X )M = F ( X ) X + Bm ( X )M

(49)

  A( X ) Bm ( X ) = A0 + θ X + B m0 + θ M θ θ

6×6

where σ ∈ R is a positive definite matrix which can be estimated and used as a tunable design parameter. Then d ( X ) ≤ g max ( X ) is used to obtain 2 g max (X)

T

= (C X ) σ (C X )

(50)

where the constant matrix C is given by



C=

I 6×6 06×6 06×6 06×6



(51)

  +  R˜ B ( X )2 g 2

max ( X )

 2 =  R˜ B + ( X ) X T C T σ C X





X  = f ( X ) + B ( X )U ( X ) + I − B ( X ) B + ( X ) v ( X )

= f ( X ) + Bm( X )M ( X ) where B m ( X ) = [ B ( X )

I 12×12 − B ( X ) B + ( X ) ], M ( X ) =

(53)  U (X)  v(X)

.

Through Eq. (49), the cost function (36) can be rewritten in a quadratic form

J=

∞ 

2

⎡ ⎢a ⎢ 11 ⎢a ⎢ 21 ⎢ F ( X ) = ⎢ a31 ⎢ ⎢ a41 ⎢ ⎣ a51

T

T



e Q e + M R M dτ

(54)

ρ 2 C T σ C + Q˜ and R =

 2 R˜ T R˜

 2

ρ I 6×6

a13 a23 a33 a43 a53 0

0 0 0 0 0 0 0 0 0 a54 0 0

0 0 0 0 0 a64

0 0 −2 0 0 0

I 6× 6 0 2 0 0 0 0

⎤ 0 0 0 0 0 0 0 a44 0 a55 0 a65

0 0 0 a45 a56 a66

0 0 0 a46 0 0

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

The detailed expression of F ( X ) is given in Appendix A. A 0 and B 0 are constant matrices. ( A 0 , B m0 ) is a controllable pair and [ F ( X ), B m ( X )] is point-wise controllable. Here, we choose A 0 = F ( X 0 ), B m0 = B m ( X 0 ) for more information about the dynamics, where X 0 = X (0) for short. Then, the new co-state equations are obtained as

∂H = F ( X ) X − B m R −1 B m ( X ) T λ ∂λ   ∞  ∂H  i λ =− =− Q + D i θ e − F T ( X )λ ∂X X =

are positive

semi-definite and positive definite weighting matrices, respectively.

(61) (62)

i =0

Once we set up A 0 , A ( X ), Q , and R, we can get the θ –D nonlinear suboptimal controller of the Coulomb tether satellite system as follows

0

where Q =

0 6× 6 a12 a22 a32 a42 a52 0

(60) (52)

As a result, the condition (52) can be satisfied with a proper selection of ρ , σ and Q˜ . According to Section 3.2, the optimal control law U ( X ) and v ( X ) can be found to minimize the cost function (36) subject to

1

where

0

Then,

(59)

U ∗( X ) v ∗( X )



= M ∗ ( X ) = − R −1 B m ( X ) T λ = − R −1 B m ( X )T Pˆ ( X , θ)e

(63)

JID:AESCTE

AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.8 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

8

Table 1 Simulation parameters used for reconfigurations. Parameter

Value

Units

m1 m2 l

150 150 3.89 × 105 2.662 × 10−6 1.991 × 10−7 15 30 1 2 6.531 × 104

kg kg km rad/s rad/s deg deg m2 m2 km

ω ωB γ δ A1 A2 X L2

The detailed process of the derivation of the θ − D controller is given in Appendix B. D 1 , D 2 , and D 3 in Eqs. (B.9)–(B.11) are chosen as the perturbation matrices. Note that only U ∗ ( X ) is the Coulomb satellite system robust suboptimal control law and v ∗ ( X ) is just an ancillary control. Though the input U ∗ ( X ) is a suboptimal controller due to the instrumental variables D and ∞θ , it isi extremely close to the optimal controller since the term i =0 D i θ is extremely small, and the variables can be used to modulate system transient responses. The θ –D method is particularly useful for solving the nonlinear optimal control problem since the suboptimal control law can be obtained in a closed form by virtue of the θ –D algorithm, which facilitates the real-time implementation since it does not demand intensive computation load.

Fig. 4. Simulation results of solar drag forces at Earth–Moon L 2 for nominal initial conditions.

4. Numerical simulations In this section, we present numerical simulations to illustrate the effectiveness of the proposed suboptimal robust controller. X L2 is defined as the X Position of L 2 point in the frame OXYZ. The paper is based on Inampudi’s work [11,12]. They had discussed the rationality of the data used for Coulomb tethered satellites, so the Coulomb tethered satellite data in the case study came from Refs. [11,12]. However, they did not consider the system orbit motion. We also use some orbit information of Earth–Moon L 2 point from CRTBP in the simulations, which has been used in many relevant papers. Table 1 provides the simulation parameters and their values. The dynamics of Coulomb satellite system is simulated using the numerical differential equation solver ODE45 in Matlab for several cases. The relative tolerance and the absolute tolerance are 10−8 and 10−10 , respectively. In this simulation, the center of mass C of the Coulomb system will be controlled to stay at the Earth–Moon L 2 point for an observation mission, and the inter-craft distance increases from 25 m to 35 m. It is assumed that the initial position errors in the x, y, and z direction are −190 km, +155.6 km, and −155.6 km, respectively. The initial libration angles are α (0) = 0.1 rad, β(0) = 0.1 rad. The nondimensional initial position and libration angular rates are x (0) = 0.1, y  (0) = −0.1, z (0) = 0, α  (0) = 0.1, β  (0) = −0.1. The desired non-dimensional position and attitude values and the final rates are set to

xf =

X L2

,

l y f = z f = α f = β f = xf = y f = zf = α f = β f = 0

(64)

Fig. 4 shows the time histories of differential solar drag force on a two-satellite Coulomb tether system at Earth–Moon L 2 point, where the satellite cross-sectional areas for nominal initial conditions are fixed. It clearly shows that solar drag force cannot be neglected in the design of the control law compared with the Coulomb force magnitude.

Fig. 5. The motion of the Coulomb Tether System (CTS) with control.

In the simulation, the actual parameter uncertainty on the cross section areas of the satellites are assumed to be

 A 1 = 0.5 A 1 ,

 A 2 = 0.3 A 2

(65)

The parameters of the uncertainty bounds are chosen to be ρ = 200, σ = 10I 6×6 . ˜ and Q in the According to the condition (40), the matrices R ˜ = diag([10, 10, 10, 10, 10, cost function (54) are selected to be: R 10]), Q = diag([2 × 106 , 106 , 106 , 105 , 105 , 109 , 6 × 103 , 6 × 103 , 6 × 103 , 6 × 103 , 6 × 103 , 6 × 108 ]). The design parameters in the perturbation matrices D i are chosen to be ki = 1, and c i = 10 (i = 1, 2). Besides, the Coulomb propulsion and electric propulsion, i.e., T X , T Y , T Z , T θ , T ϕ , T L can be obtained as follows

T X = mω2lU 1 ,

T Y = mω2lU 2 ,

T α = bmω2 U 4 ,

T β = bmω2 U 5 ,

T Z = mω2lU 3 , T L = bmω2lU 6

The charges q1 , q2 in magnitude across the satellites are obtained by Eqs. (25)–(27). To test the robustness of the controller, the external disturbance due to the differential solar radiation pressure effects is added to the system simulation. The total simulation time is about two orbital periods (54.6372 days). Figs. 5 and 6 reflect the orbit motion of the Coulomb tether satellite system during the reconfiguration. Fig. 5 shows the actual trajectory of the barycenter of the system with the proposed controller in about two orbital periods. It is observed that the system

JID:AESCTE AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.9 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

9

Fig. 6. The time history of the CTS position error. Fig. 9. The time history of the control input for the X , Y , Z axes.

Fig. 7. The time history of the libration angles

α and β . Fig. 10. The time history of the control input for the libration angles.

Fig. 8. The time history of the distance between the satellites.

reaches the desired position and is kept at L 2 point with small position errors under the external disturbances. As can be seen from the time history of the orbit position error in Fig. 6, after about 30 days, the position error goes to zero, which indicates that the center of mass of the system moves to L 2 point according to the requirement with a good station-keeping performance. Fig. 7 shows the Coulomb satellite attitude motion during the Coulomb satellite system reconfiguration. The in-plane pitch angle α and out-of-plane motion angle β asymptomatically converge to zero in less than one orbit period (27.3186 days). It is seen from Fig. 8 that the separation distance L d between the satellites expands from an initial 25 m to a final 35 m. The reconfiguration is also finished in one orbit. These results demonstrate the precise and robust performance of the proposed controller under the external disturbance due to the differential solar radiation pressure effects.

Fig. 11. The time history of the control input for the length L d .

Figs. 9 and 10 present the control forces for the motion of position of the system and attitude libration angles. The initial control efforts are relatively large in order to drive the system to the desired position and attitude quickly. They decrease rapidly after the desired position and attitude are achieved. Very small forces are required to maintain the position and configuration of the system due to the disturbance. As Fig. 9 shows, the thrusts along the X , Y , and Z directions drop quickly on the first 5 days due to the initial position errors and eventually remain near zero. The amplitude of thrusts required by the controller is less than 1.5 mN, which is reasonable and can be implemented through current low-thrust engines. Figs. 10 and 11 show inertial thruster usage for in-plane

JID:AESCTE

AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.10 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

10

Council (CSC) for financially supporting them as visiting scholars at Concordia University. The authors would also like to thank the associate editor and reviewers for their very constructive comments and suggestions which have helped greatly improve the quality and presentation of the paper. Appendix A. Terms in Eq. (60)

a11 =

1−ν

+ −

Fig. 12. The time history of the satellite charge q1 .

and out-of-plane control and the Coulomb force for longitudinal control. It can be observed that Coulomb control and transverse control (micro-thrusters) forces are on the order of μN. Transverse control can be implemented either by using Colloid or PPT microthrusters. The charge on the second satellite will be equal and opposite to that of satellite 1. Fig. 12 shows the time histories of satellite control charge q1 for one of the satellites. At the end of the reconfiguration, the Coulomb control charge should explicitly approximate the equilibrium value. The maximum of charge level is about 2.3 μC, which is easily implementable by the charge emission devices. From these simulation results, we can obtain that the developed suboptimal robust control scheme based on the θ –D method is valid for the Coulomb satellite system control problem. It not only makes the two-satellite Coulomb formation reconfigure successfully, but also keeps the system staying steadily at Earth– Moon L 2 point with relatively low energy consumption. 5. Conclusion This paper investigates the dynamics and reconfiguration control problem of a two-satellite Coulomb tether formation in the presence of differential solar drag near Earth–Moon libration point. An indirect robust control scheme is designed via an optimal control formulation. This indirect robust control addressed the coupling motion of attitude and orbit in one unified optimal control framework. The θ –D technique is applied to solve the resultant nonlinear control problem. Numerical simulation results demonstrate the capacity of the proposed controller in terms of precise and robust performance under external disturbance. This investigation restricts attention to a two-body Coulomb satellite system reconfiguration and station-keeping problem at libration points. For such study, the advantage of using coupling motion of attitude and orbit is highlighted. It is beneficial to future works. The reconfiguration and station-keeping problem presented here can be extended to periodic orbits of the Coulomb satellite system in the CRTBP. Some similar station-keeping strategies for the Coulomb satellite system moving along the periodic orbits will be obtained. Also, the proposed controller can be extended to reconfigure three-craft formations.

+ a12 =

a21 =

The authors declare that they do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted. Acknowledgements This present work was supported by the National Natural Science Foundation of China (NSFC) under Grant 61304005. Jing Huang and Gang Liu would like to acknowledge China Scholarship

ν r23

(1 − ν )εb(1 + x) S 2 (cos χ1 ) xr15 2

3νεb cos β cos2 α r25

νεb S 2 (cos χ2 )

−

(A.1)

r25

3νεb cos2 β sin α cos α

(A.2)

r25 r15

+

3νεb sin β cos β cos α r25

(A.3)

3(1 − ν )εb(1 + x) cos2 β sin α cos α xr15

−

3νεb cos2 β sin α cos α

(A.4)

r25 1−ν

ν

−

r13

r23

+

3(1 − ν )εb cos2 β sin2 α

(1 − ν )εb S 2 (cos χ1 ) r15

r15

+

3νεb cos2 β sin2 α r25

νεb S 2 (cos χ2 )

(A.5)

r25

3(1 − ν )εb sin β cos β sin α r15

+

3νεb sin β cos β sin α r25

(A.6)

3(1 − ν )εb(1 + x) sin β cos β cos α xr15

+

3νεb sin β cos β cos α

(A.7)

r25

3(1 − ν )εb sin β cos β sin α r15

a33 = −

−

Conflict of interest statement

−

xr15

3(1 − ν )εb sin β cos β cos α

−

a32 =

xr13

3(1 − ν )εb(1 + x) cos2 β cos2 α

a22 = 1 −

a31 =

1−ν

r15

+

a23 =

−

r13

3(1 − ν )εb cos2 β sin α cos α

+ a13 =

1−ν

+1−

x

1−ν r13

−

ν r23

+

+

3νεb sin β cos β sin α r25

(A.8)

3(1 − ν )εb sin2 β

(1 − ν )εb S 2 (cos χ1 ) r15

r15

+

3νεb sin2 β r25

−

νεb S 2 (cos χ2 ) r25 (A.9)

a41 = − a42 =

3(1 − ν )(1 + x)2 sin α cos α

3(1 − ν ) 

xr15

r15

+

3ν  r25

−

3ν x sin α cos α r25

(A.10)

(1 + x) cos2 α − (1 + x) sin2 α + y sin α cos α

x cos2 α − x sin2 α + y sin α cos α





(A.11)

JID:AESCTE AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.11 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

a43 =

3(1 − ν )  r15

+

3ν r25

−(1 + x) sin α tan β + y cos α tan β

(−x sin α tan β + y cos α tan β)

a44 = 2β  tan β +

2ε 

a52 =

(A.14) (A.15)

ε

−

3(1 − ν )(1 + x)2 sin β cos β cos2 α 3ν x sin β cos β cos2 α

(A.16)

r25 r15

−2(1 + x) sin β cos β sin α cos α

r2

− y sin2 β sin α + z sin β cos β



sin β cos β

T 1 A 0 − B m0 R





a66 = ε β

=0

(B.4)

T B m0 T0 T



+

A ( X )T 0

θ



A 0T

− T 0 B m0 R

+ T 0 B m0 R

−1

T −1 B m ( X )





(B.5) T1

T0

θ

B m ( X ) −1 T R B m0 T 0 − D 1

θ

T B m0

(B.6)





+

2ν P 2 (cos χ2 ) r23

T 1 A( X )

2

+ cos β

ε cos2 β



θ θ

T1

θ

θ

−1

T B m0 T1

+ T 1 B m0 R

θ

T −1 B m ( X )

θ

T0

B m ( X ) −1 T R B m0 T 0 − D 2

(B.7)

θ

+

T n −1 A ( X )

θ n −2 

T

j

j =0

+

n −1 





A T ( X ) T n −1

θ

T B m ( X ) −1 B m (X) R T n −2 − j

θ

T

−

j

θ

B m0 R −1

T Bm (X)

θ

j =0

+

(A.22)

n −1 



+

B m ( X ) −1 T R B m0 T n−1− j

θ

T T j B m0 R −1 B m0 T n− j − D n

(B.8)

j =1

(A.23)

The perturbation matrix D i (i = 1, 2, . . . , n) is constructed as follows:



D 1 = k1 e −c1 t −

T 0 A( X )

θ

+ T 0 B m0 R (B.1)

D 2 = k2 e −c2 t −

i =0

+ T0

From the optimal control theory [16], for an infinite-time optimal problem, it can be inferred that t →∞

T Bm (X)

+ T 0 B m0 R −1

T B m ( X ) −1 T B m ( X ) −1 B m (X) R B m0 T 1 + T 0 R T0



P ( X , θ) ≈ lim P ( X , θ) = const

θ



=−

(A.24)

T i ( X , θ)θ i X = P ( X , θ) X

A T ( X )T 1

−

T T T n A 0 − B m0 R −1 B m0 T 0 + A 0T − T 0 B m0 R −1 B m0 Tn

Assuming a power series expansion of λ ∞ 

Diθ

.. . 

Appendix B. Derivation of the θ –D controller [29]

λ=

+ T0

(A.21)



−

θ

+ T1

(A.20)

ε

a65 = 2 + α 

T 0 A( X )

−1

+ T 1 B m0 R

(A.19)

2ε 

r13

 i

T ( X ) P ( X , θ) P ( X , θ) B m ( X ) R −1 B m

i =0



(A.18)

a55 = −2 sin β cos β − α  sin β cos β

a64 =

∞ 

2

T T 0 A 0 + A 0T T 0 − T 0 B m0 R −1 B m0 T0 + Q = 0

(A.17)

 − (1 + x) sin2 β cos α − y sin2 β sin α + z sin β cos β 3ν  + 5 x cos2 β cos α + y cos2 β sin α − x sin2 β cos α

2(1 − ν ) P 2 (cos χ1 )

Q +

+ T0

r1

a56 = −



1

P ( X , θ) −

θ

In order to determine P (x, θ), let the coefficients of powers of

=−

 − y sin β cos β sin2 α 3(1 − ν )  (1 + x) cos2 β cos α + y cos2 β sin α a53 = 5

β

T

θ i in Eq. (B.4) be zero:

r2

a54 = −

2

A( X )

T T T 2 A 0 − B 0 R −1 B m0 T + A 0T − T 0 B m0 R −1 B m0 T2

 − y sin β cos β sin2 α 3ν  + 5 −2x sin β cos β sin α cos α

1

1

=−

xr15

3(1 − ν ) 

+

(A.13)

ε

2

a51 = −

A0 + θ (A.12)

a45 = 2 tan β a46 =





11

+ T1

Therefore, the derivative of Eq. (B.1) is given by

  T λ = P ( X , θ) X  = P ( X , θ) F ( X ) X − B m ( X ) R −1 B m ( X )λ (B.3) The algebra Riccati equation and the co-state vector equation can be derived from Eqs. (56), (B.1), and (B.3) as

θ θ

A T ( X )T 0

θ

T0 + T0

−

B m ( X ) −1 T R B m0 T 0

 (B.9)

θ

A T ( X )T 1

θ

T T1 + T 1 B m0 R −1 B m0

T B m ( X ) −1 T B m ( X ) −1 B m (X) R B m0 T 0 + T 0 R T0

θ

θ

−1

T B m0 T1

+ T 1 B0 R 

θ

T −1 B m ( X )

θ

B ( X ) −1 T R B m0 T 0

.. . D n = kn e

T −1 B m ( X )

T 1 A( X )

+ T 1 B m0 R

(B.2)

−

(B.10)

θ

 −cn t

−

T n −1 A ( X )

θ

T0

−

A T ( X ) T n −1

θ

JID:AESCTE

AID:3070 /FLA

[m5G; v 1.134; Prn:9/07/2014; 8:14] P.12 (1-12)

J. Huang et al. / Aerospace Science and Technology ••• (••••) •••–•••

12

+

n −2 

T

j =0

+

n −1 

T B m ( X ) −1 B m (X) R T n −2 − j

j

θ

T

j

θ

B 0 R −1

T Bm (X)

θ

j =0

+

n −1 

[7]



+

[8]

B m ( X ) −1 T R B m0 T n−1− j

θ

[9]

 T T j B m0 R −1 B m0 T n− j

(B.11)

[10]

where ki and c i > 0, i = 1, 2...n are design parameters. Then, Eqs. (B.5)–(B.8) become

[11]

j =1

−

T n −1 A ( X )

−

θ +

n −1 

A T ( X ) T n −1

θ

T

B m0 R

j

n −1 

T B m ( X ) −1 B m (X) Tj R T n −2 − j

θ

j =0

T −1 B m ( X )

θ

j =0

+

+

n −2 

θ

[13]



B m ( X ) −1 T R B m0 T n−1− j +

[14]

θ

[15]

T T j B m0 R −1 B m0 T n− j − D n

j =1



[16]

= ζi (t ) − +

n −2 

T n −1 A ( X )

θ

n −1  j =0

+

n −1 

−

A T ( X ) T n −1

[17]

θ [18]

T B m ( X ) −1 B m (X) Tj R T n −2 − j

θ

j =0

+

[12]

T

j

B m0 R −1

θ

T Bm (X)

θ

[19]



+

[20]

B ( X ) −1 T R B m0 T n−1− j

θ

 T T j B m0 R −1 B m0 T n− j

[21]

(B.12)

j =1

[22]

where ζi (t ) = 1 − ki e −c i t .ζi (t ) is chosen to satisfy the convergence and stability conditions [29]. The first three terms T 0 , T 1 , and T 2 in P ( X , θ ) are used to get the solution. Therefore,

Pˆ ( X , θ) ≈ T 0 + T 1 ( X , θ)θ + T 2 ( X , θ)θ 2

(B.13)

[23]

[24] [25]

References [1] J. Berryman, Analytical and numerical analysis of static Coulomb formations, Master’s Thesis, Aerospace and Ocean Engineering Department, Virginia Polytechnic Institute and State University, Blacksburg, VA, December 2005. [2] J. Berryman, H. Schaub, Analytical charge analysis for 2- and 3-craft Coulomb formations, AIAA J. Guidance Control Dyn. 30 (6) (Nov.–Dec. 2007) 1701–1710. [3] Z. Cai, X. Li, Z. Wu, Deployment and retrieval of a rotating triangular tethered satellite formation near libration points, Acta Astronaut. 98 (2014) 37–49. [4] Y. Chen, R. Huang, L. He, Dynamical modelling and control of space tethers: a review of space tether research, Nonlinear Dyn. (2014) 1–23. [5] J.H. Cover, W. Knauer, H.A. Maurer, Lightweight reflecting structures utilizing electrostatic inflation, US Patent 3,546,706, Oct. 1966. [6] R. Farley, D. Quinn, Tethered formation configurations: meeting the scientific

[26]

[27] [28] [29] [30]

objectives of large aperture and interferometric science, in: Proceedings of 2001 AIAA GNC Conference, 2001-4770. S. Gates, Dynamic model for a multi-tethered space-based interferometer, Report NRL/MR/ 8230-00-8475, Naval Research Laboratory, August 2000. D. Hastings, H. Garrett, Spacecraft-Environment Interactions, Cambridge University Press, 2004. H. Huang, L. Yang, Y. Zhu, Dynamics and relative equilibrium of spacecraft formation with non-contacting internal forces, Proc. Inst. Mech. Eng. G J. Aerosp. Eng. 228 (7) (June 2013) 1171–1182. R. Inampudi, Optimal reconfigurations of two-craft Coulomb formation in circular orbits, in: AAS/AIAA Astrodynamics Specialist Conference Guidance, Toronto, Canada, Aug. 2010, pp. 2–5. R. Inampudi, H. Schaub, Orbit radial dynamic analysis of two-craft Coulomb formation at libration points, in: AIAA/AAS Astrodynamics Specialist Conference, Toronto, Ontario Canada, Aug. 2–5, 2010. R. Inampudi, H. Schaub, Orbit-radial control of a two-craft Coulomb formation about circular orbits and libration points, in: 4th International Conference on Spacecraft Formation Flying Missions & Technologies St-Hubert, Quebec, May 18–20, 2011. R. Jones, Optimal reconfigurations of two-craft Coulomb formations along manifolds, Acta Astronaut. 83 (2013) 108–118. L.B. King, G.G. Parker, S. Deshmukh, J.-H. Chong, Study of inter-spacecraft Coulomb forces and implications for formation flying, AIAA J. Propuls. Power 19 (3) (May–June 2003) 497–505. L.B. King, G.G. Parker, S. Deshmukh, J.H. Chong, Spacecraft formation-flying using inter-vehicle Coulomb forces, Tech. Rep., NASA/NIAC, Jan. 2002, http://www. niac.usra.edu. D.E. Kirk, Optimal Control Theory: An Introduction, Dover Publications Inc, Mineola, 2003. F. Lin, Robust Control Design: An Optimal Control Approach, John Wiley & Sons Ltd, Hoboken, 2007. G. Liu, J. Huang, G.F. Ma, Nonlinear dynamics and station-keeping control of a rotating tethered satellite system in Halo Orbits, Chin. J. Aeronaut. 26 (5) (2013) 1227–1237. C.R. McInnes, Solar Sailing: Technology, Dynamics and Mission Applications, Springer-Praxis, Chichester, UK, 1999. A.K. Misra, J. Bellerose, V.J. Modi, Dynamics of a tethered system near the Earth–Moon Lagrangian points, in: Proceedings of the 2001 AAS/AIAA Astrodynamics Specialist Conference, Quebec City, Canada, in: Advances in the Astronautical Science, vol. 109, 2002, pp. 415–435. N. Murdoch, D. Izzo, C. Bombardelli, et al., The electrostatic tractor for asteroid deflection, in: 58th International Astronautical Congress, 2008, Paper IAC-08-A3.I.5. A. Natarajan, A study of dynamics and stability of two-craft Coulomb tether formations, Ph.D. Dissertation, Aerospace and Ocean Engineering Department, Virginia Polytechnic Institute and State University, Blacksburg, VA, May 2007. A. Natarajan, H. Schaub, Hybrid control of orbit-normal and along-track two-craft Coulomb tethers, in: AAS/AIAA Space Flight Mechanics Meeting, No. 07-193, Jan.–Feb. 2007. J. Peláez, M. Sanjurjo, Generator regime of self-balanced electrodynamic bare tethers, J. Spacecr. Rockets 43 (6) (2006) 1359–1369. L. Pettazzi, D. Izzo, Self-assembly of large structures in space using intersatellite Coulomb forces, in: 56th International Astronautical Congress, Fukuoka, Japan, Oct. 17–21, 2005, Paper IAC-06-C3.4/D3.4.07. M. Sanjurjo-Rive, F.R. Lucas, J. Peláez, On the dynamics of a tethered system near the collinear libration points, in: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, Hawaii, Aug. 18-21, 2008, pp. 18–21, AIAA 2008-7380. H. Vasavada, H. Schaub, Analytic solutions for equal mass 4-craft static Coulomb formation, J. Astronaut. Sci. 56 (1) (Jan.–Mar. 2008) 7–40. B. Wong, A.K. Misra, Dynamics of Lagrangian point multi-tethered satellite systems, J. Astronaut. Sci. 53 (3) (2005) 221–250. M. Xin, S.N. Balakrishnan, A new method for suboptimal control of a class of non-linear systems, Optim. Control Appl. Methods 26 (3) (2005) 55–83. M. Xin, H. Pan, Indirect robust control of spacecraft via optimal control solution, IEEE Trans. Aerosp. Electron. Syst. 48 (2012) 1798–1809.