arXiv:1603.05244v1 [math.OC] 16 Mar 2016
Three-Dimensional Multi-Tethered Satellite Formation
with the Elements Moving Along Lissajous Curves
D. Yarotsky∗
V. Sidorenko†
D. Pritykin‡
Abstract
This note presents a novel approach to maintain three-dimensional multi-tethered
satellite formation in space. For a formation consisting of a main body connected by
tethers with several deputy satellites (the so-called “hub-and-spoke” configuration) we
demonstrate that under proper choice of the system’s parameters the deputy satellites
can move along Lissajous curves in the plane normal to the local vertical with all
tethers stretched; the total force due to the tension forces acting on the main satellite
is balanced in a way allowing it to be in relative equilibrium strictly below or strictly
above the system’s center of mass. We analyze relations between the system’s essential
parameters and obtain conditions under which the proposed motion does take place.
We also study analytically the motion stability for different configurations and whether
the deputy satellites can collide or the tethers can entangle. Our theoretical findings
are corroborated and validated by numerical experiments.
Keywords: Tethered Satellite System, Satellite Formation, Dynamics, Control,
Stability
1
Introduction
Three-dimensional satellite formations are often discussed in connection with multi-point
measurements needed for atmospheric, geodetic or plasma physics studies. To simplify control strategies and to minimize fuel consumption, tethers can be used to maintain desired
relative positions of satellites in the formation flying. For the first time three-dimensional
multi-tethered formations were discussed probably by [6], who proposed double-pyramid configurations. It seems that the most straightforward way to obtain a multi-tethered formation
is to deploy from the main satellite several tethers with deputy satellites at their ends. To
∗
Institute for Information Transmission Problems, Russian Academy of Sciences, Bolshoy Karetny per.
19, 127051 Moscow, Russian Federation
†
Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, Miusskaya Sq., 4, 125047
Moscow, Russian Federation
‡
Moscow Institute of Physics and Technology, Institutskiy per., 9, 141700, Dolgoprudny, Moscow Region,
Russian Federation
1
specify such formations [13] introduced the term “hub-and-spoke”. Behavior of “hub-andspoke” multi-tethered formations has been studied for different dynamical environments: in
circular orbit [1, 4], in elliptic orbit [5], in halo-orbit [7, 18] and near collinear Lagrangian
points [17].
To keep the tethers taut the combination of rotation with gravity-gradient forces is usually
proposed. Among other opportunities the relatively new concept of the Tethered Coulomb
Structure (TCS) is worth mentioning [12,14,15]. In this case the satellites are electrostatically
charged to produce repulsive forces between them. Nevertheless, it looks as if Coulomb
repulsive forces can be used to prevent the slack of short enough tether: in [12, 14, 15] the
discussed length is 10 m by order of magnitude.
To give an idea of our approach we begin with the system of two bodies connected by
a single tether aligned along the local vertical; the mass center of the system moves in a
circular orbit. Small oscillations of this system around local vertical are a combination of
in-plane and out-of-plane natural oscillations with incommensurable frequencies [8, 16]. The
motion of end bodies in these oscillations can be roughly described as a motion along curves
densely filling certain areas on planes normal to the local vertical. Then let us consider the
degenerate “hub-and-spoke” configuration in the relative equilibrium with all tethers aligned
along the vertical (Fig. 1, left) The displacement of a single deputy satellite from the relative
equilibrium position causes oscillations whose frequencies differ from those inherent in the
preceding case. With the proper tuning of the system’s parameters these frequencies can be
made commensurable resulting in the motion of the deputy satellite along a Lissajous curve
in the plane, normal to the local vertical (Fig. 1, right). Naturally all other deputy satellites
can also be put in motion along similar curves. The choice of initial conditions allows to
avoid collisions among them and to ensure the balance of the tension forces applied to the
main satellite so as to preserve its relative equilibrium.
We suppose that the described structure can be useful for some applications or at least
become the starting point for the development of new approaches to maintain 3D multitethered satellite formations in space.
In Section 2 we start our study with a simplified dynamical model of multi-tethered
formation (point masses + weightless tethers). In Section 3 we consider the linearized dynamics near vertical equilibrium and describe the oscillations in the system. In Section 4 we
consider deputy satellite formations moving along Lissajous curves. In Section 5 we present
the results obtained by numerical simulation of the system’s dynamics.
2
Deputy satellite dynamics at small deviations from
the relative equilibrium
As mentioned in the Introduction, we will consider the multi-tethered satellite formation
comprising N +1 bodies: the main satellite C of mass mC and N deputy satellites D1 , ..., DN
(each of mass mD ) linked to the main satellite (but not to each other) by identical extensible
tethers; the tethers’ masses are ignored.
2
D2
D2
D3
D1
D3
D1
O y
x
z
C
C
Figure 1: Multi-tethered formation consisting of the main body C and three deputy
satellites D1 , D2 , D3 . On the left the system is in relative equilibrium with all tethers
aligned along the local vertical. On the right the deputy satellites move along a Lissajous curve (see also the animation provided as Electronic Supplementary Material).
To write down the equations of motion, we shall introduce a Local Vertical Local Horizontal (LVLH) reference frame Oxyz, centered on the position of the system’s center of mass
(CoM) in its nominal orbital motion: Oz axis is aligned with the local vertical and oriented
towards Earth’s center, Ox runs tangentially to the orbit in the direction of the CoM motion,
and Oy axis is directed along the normal line to the orbit plane (Fig. 1). It is assumed that
nominally the system’s center of mass moves in circular orbit with the mean motion ω0 . The
unit vectors in the directions of the axes x, y, z will be denoted as ex , ey , ez respectively.
In the LVLH frame the motion of any satellite in the considered formation can be described by the Hill-Clohessy-Wiltshire (HCW) equations [9]:
ẍ = 2ω0 ż + Tx /m,
ÿ = ω02 y + Ty /m,
z̈ = −2ω0 ẋ + 3ω02 z + Tz /m,
(1)
where Tx , Ty , Tz denote the components of the (total) tether tension T applied to a given
body of mass m.
We shall adopt the usual visco-elastic model of massless tether, hence the tension force
applied to the ith deputy satellite is
d
(rC − ri )
Ti = 1(|rC −ri |>l0 ) · k(|rC − ri | − l0 ) + b |rC − ri |
,
(2)
dt
|rC − ri |
3
where l0 is the slack tether length, k the elastic coefficient, b the damping coefficient, and
(
1, if |rC − ri | > l0
1(|rC −ri |>l0 ) =
0, otherwise
The vertical equilibrium of the system takes place if all tethers rest (ṙC = ṙi = 0) aligned
along the local vertical:
yi∗ = yC∗ = 0,
−1
mC ∗
N m D + mC
mD ω02
∗
zi = −
.
z = l0
−3
N mD C
mC
k
x∗i = x∗C ,
(3)
Here and below equilibrium quantities are marked by the asterisk *. In contrast to the
y and z coordinates, the x coordinate in this configuration, though shared by all satellites,
is arbitrary, since the equations are invariant with respect to translations along the orbit
(i.e., along the x axis in the approximation of the orbital dynamics, provided by the HCW
equations).
For the above configuration to be valid the denominator in the last formula must be
positive, i.e. the tethers must be sufficiently rigid to counteract the microgravity:
mC mD
.
(4)
N m D + mC
Denoting the tether’s length in the equilibrium configuration by l∗ = |r∗C − r∗i |, we obtain
a useful relation
k > 3ω02
l∗ − l0
3ω 2 mC mD
= 0
.
l∗
k N m D + mC
(5)
Denote
3ω 2
mC mD
l∗ − l0
= 0 m r , mr =
.
l∗
k
N m D + mC
The tension forces in all the tethers in the equilibrium configuration (3) have the same
value T∗ given by the obvious formula
λ∗ =
T∗ = k(l∗ − l0 )ez .
(6)
It is not difficult to derive an approximate expression for the tension forces in case of small
displacements of the main satellite and the ith deputy satellite with respect to equilibrium
(3):
Ti =T∗ + kλ∗ (∆xC − ∆xi )ex + kλ∗ (∆yC − ∆yi )ey +
+ [k(∆zC − ∆zi ) + b(∆żC − ∆żi )] ez .
Here ∆xC , ∆yC , ∆zC are the components of the main satellite’s displacements, and ∆xi ,
∆yi , ∆zi are the components of the ith deputy satellite displacements.
4
With the linearized tension the HCW equations read
k
(∆xC − ∆xi ),
mD
k
(∆yC − ∆yi ),
∆ÿi = −ω02 ∆yi + λ∗
mD
b
k
(∆zC − ∆zi ) +
(∆żC − ∆żi )
∆z¨i = −2ω0 ∆ẋi + 3ω02 ∆zi +
mD
mD
∆ẍi = 2ω0 ∆z˙i + λ∗
(7)
for the deputy satellites, and
N
k X
∆x¨C = 2ω0 ∆z˙C − λ∗
(∆xC − ∆xi ),
mC i=1
∆y¨C =
−ω02 ∆yC
N
k X
− λ∗
(∆yC − ∆yi ),
mC i=1
∆z¨C = −2ω0 ∆x˙C +
3ω02 ∆zC
(8)
N
N
b X
k X
(∆zC − ∆zi ) −
(∆żC − ∆żi )
−
mC i=1
mC i=1
for the main satellite.
3
Decoupling the equations of motion
The system of equations derived above can be split up into three independent groups by
taking appropriate linear combinations, separately for x, y and z components:
1. By taking the sum of the equation for the main satellite with weight mC and all the
respective equations for deputy satellites with weights mD we obtain a triple of scalar
equations for
!
N
X
1
mC ∆rC + mD
∆ri ,
N m D + mC
i=1
i.e., for the motion of the whole system’s center of mass.
2. By taking the sum of the equations for deputy satellites with coefficients 1/N and
subtracting from it the equation for the main satellite we obtain a triple of scalar
equations for
N
1 X
∆ri − ∆rC ,
N i=1
i.e. describing the relative motion between the main satellite and the auxiliary satellites’ center of mass.
5
3. Finally, if we take any linear P
combination of equations for deputy satellites with
N
some coefficients σi such that
i=1 σi = 0, we obtain a triple of scalar equations
PN
for i=1 σi ∆ri , which can be interpreted as a partial description of the relative motion between deputy satellites. To obtain the full description, we need to consider all
N − 1 linearly independent such assignments of coefficients, thus giving the total of
3(N − 1) scalar equations. The simplest example of the set with suitable combinations
of coefficients is
{−1, 1, 0, . . . , 0}, {0, −1, 1, 0, . . . , 0}, . . . , {0, . . . , −1, 1}.
Physically it means that we use relative displacements ∆ri+1 − ∆ri (i = 1, N − 1) of
consecutive objects to describe the dynamics of the subsystem composed of the deputy
satellites. One more opportunity is provided by the set
1
1
1
1
1
1
,..., − ,...,1 − ,−
,
1 − ,− ,...,−
N
N
N
N
N
N
characterizing relative displacements of the deputy satellites D1 , . . . , D(N −1) with respect to the center of mass of all deputy satellites.
Clearly, the collection of these three groups of equations is equivalent to the original
system of 3(N + 1) scalar equations for the main and deputy satellites. Let us deal with the
three groups one by one.
3.1
System CoM motion equations
The motion of the system’s center of mass is described by free HCW equations [9]. It is
well-known that these equations describe oscillations with frequency ω0 in the y component,
and oscillations with frequency ω0 combined with a linear drift in the orbital xz plane.
3.2
Relative motion between the main satellite and the deputy
satellites’ center of mass
Let ∆x, ∆y, ∆z denote the scalar components of
1
N
PN
i=1
∆ri − ∆rC . Then using relation (5)
k
∆x = 2ω0 ∆ż − 3ω02 ∆x,
mr
k
∆ÿ = −ω02 ∆y − λ∗ ∆y = −4ω02 ∆y,
mr
b
k
∆z −
∆ż.
∆z̈ = −2ω0 ∆ẋ + 3ω02 ∆z −
mr
mr
∆ẍ = 2ω0 ∆ż − λ∗
(9)
The equation for ∆y is independent from the equations for ∆x, ∆z and describes harmonic
oscillations with the frequency 2ω0 . To analyze the remaining equations for ∆x, ∆z we
consider the corresponding first-order system
6
0
1
0
∆x
2
d ∆ẋ −3ω0
0
0
=
0
0
0
∆z
dt
2
0
−2ω0 3ω0 −
∆ż
k
mr
0
∆x
2ω0 ∆ẋ
1 ∆z
− mbr
∆ż
.
(10)
Let us denote by A the matrix on the right-hand side of the system (10). To examine
the stability property of this system, we write down the characteristic equation
b 4
3bω02
k
k
4
2
2
2
2
det(A − ρI4 ) = ρ +
ρ +
+ 4ω0 ρ +
ρ + 3ω0
− 3ω0 = 0.
mr
mr
mr
mr
The symbol I4 is used here for the identity matrix of the fourth order. Applying the
Routh-Hurwitz stability criterion [10], it is not difficult to establish that all eigenvalues ρ of
A belong to the left half-plane Reρ < 0 iff the condition (4) is satisfied and b > 0.
Summarizing, under the stability assumption (4) the center of mass of deputy satellites
oscillates with respect to the main satellite; in the linear approximation the energy dissipation
in tethers does not affect the oscillations along the y axis.
3.3
Relative motion between deputy satellites
Let ∆x, ∆y, ∆z denote the scalar components of
P
ficients σi subject to N
i=1 σi = 0. Then
PN
i=1
σi ∆ri with some assignment of coef-
k
∆x,
mD
k
∆ÿ = −ω02 ∆y − λ∗
∆y,
mD
∆ẍ = 2ω0 ∆ż − λ∗
∆z̈ = −2ω0 ∆ẋ + 3ω02 ∆z −
(11)
k
b
∆z −
∆ż.
mD
mD
The equation for the y component describes oscillations with frequency
r
r
4mC + N mD
k
2
= ω0
.
ωy = ω0 + λ∗
mD
mC + N m D
The remaining equations for ∆x, ∆z can be written in the first-order form as
0
1
0
0
∆x
∆x
k
0
0
2ω0
d
∆ẋ .
∆ẋ = −λ∗ mD
0
0
0
1 ∆z
dt ∆z
k
2
∆ż
∆ż
0
−2ω0 3ω0 − mD − mbD
The matrix A1 on the right-hand side of (13) has the characteristic polynomial
b 3
k
λ∗ kb
λ∗ k
(λ∗ + 1)k
4
2
2
2
det(A1 − ρI4 ) = ρ +
ρ +
+ ω0 ρ + 2 ρ +
− 3ω0 .
mD
mD
mD
mD mD
7
(12)
(13)
Applying again Routh-Hurwitz stability criterion, we establish that the system (13) is
asymptotically stable iff b > 0 and
k ≥ 3ω02 mD .
(14)
If the condition (14) is fulfilled, but the dissipation is absent (b = 0), then the spectrum
of A1 is purely imaginary. If the condition is violated, matrix A1 has a positive eigenvalue.
Note that the stability condition (14) is stronger than the earlier condition (4).
Assuming b = 0, in the limit of large rigidity k the linearized dynamics of ∆x, ∆z
approximately decouples into independent oscillations of ∆x and ∆z with frequencies
r
mD ω02
λ∗ k
−2
1−2
=
ωx =
+O k
mD
k
r
3mC
mD ω02
−2
(15)
1−2
+O k
=
ω0 ,
N m D + mC
k
r
k
mD ω02
−2
.
1+
+O k
ωz =
mD
2k
3.4
Lyapunov function
Our solution of the linearized dynamics implies, in particular, that in the linear approximation of tether tension the vertical equilibrium is stable. This conclusion can also be shown in
a stronger sense – without linearization of tether tension – by directly providing a Lyapunov
function. Specifically, let r′ = r − rCoM , ṙ′ = ṙ − ṙCoM denote position and velocity relative
to the system’s center of mass. Consider the energy E = T + V of relative motion, where
N
mD X ′ 2 mC ′ 2
ṙ
ṙ +
T =
2 i=1 i
2 C
and
N
2
mD ω02 X ′ 2
mC ω02 ′ 2
2
V =
(y i − 3z ′ i ) +
y C − 3z ′ C +
2 i=1
2
N
+
kX
2
1(|r′ C −r′ i |>l0 ) · (|r′ C − r′ i | − l0 ) .
2 i=1
If b = 0, then Ė = 0, otherwise Ė ≤ 0. Taylor expansion of V near the equilibrium yields
8
N
m ω2
mD ω0 2 X ′ 2
C 0
′ 2
′ 2
′2
∆yC − 3∆zC +
∆yi − 3∆zi +
V = Vequilibrium +
2
2
i=1
N
io
h
k Xn
2
2
2
+
(∆zi′ − ∆zC′ ) + λ∗ (∆x′ i − ∆x′ C ) + (∆x′ i − ∆x′ C )
2 i=1
(16)
3
+ O(|∆r′ | )
Taking into account the identity
′
∆r C
N
mD X ′
∆r i ,
=−
mC i=1
one can check that this form is positive definite exactly if stability conditions (4) and (14)
hold.
4
Motion of deputy satellites along Lissajous curves
Having described oscillations of the “hub-and-spoke” system near the equilibrium, we consider now the possibility of the satellites moving in such a way that the system elements
(satellites and tethers) never collide. Given our assumption of small deviations from the
vertical equilibrium, we formulate this as the requirement that the projections (xi , yi ) of
deputy satellites to the xy plane never come close to each other.
Results of Section 3.3 imply that the position of a deputy satellite relative to the center
of mass of all deputy satellites,
N
1 X
rk ,
r′ i = ri −
N k=1
oscillates with frequency
r
4mC + N mD
mC + N m D
in the y direction and, for sufficiently rigid tethers, with frequency
r
3mC
ωx ≈ ω0
N m D + mC
ωy = ω0
(17)
(18)
in the x direction (with a more accurate value
P given by (15)). Since, as shown in Section
3.2, the center of mass of deputy satellites N1 N
k=1 rk performs independent oscillations, we
may ignore these latter and assume without loss of generality that
N
N
1 X
1 X
xk =
yk ≡ 0
N k=1
N k=1
9
(19)
at all times.
By linearity, for each pair i, j of deputy satellites their relative position ri − rj also
oscillates with the same frequencies ωx , ωy . If the frequencies are incommensurate, i.e.,
ωx /ωy is irrational, then the trajectory of the oscillation is aperiodic and it comes arbitrarily
close to the origin, i.e. the two satellites come arbitrarily close to each other. We are thus
naturally led to consider commensurate oscillations:
ωx
p
= ,
ωy
q
(20)
where p, q are co-prime natural numbers. In this case the two oscillations have a common
period of
2πp
2πq
TL =
=
.
ωx
ωy
It is convenient to introduce the non-dimensional time
τ=
t
.
TL
The xy-motion of a single deputy satellite can then be written as
x = x0 sin (2πpτ + ϕx ) , y = y0 sin(2πqτ + ϕy ),
(21)
with some initial phases ϕx , ϕy , and has a period 1. The trajectory of this motion is known
as a Lissajous curve [3, Sec. 25]
Substituting expressions (17), (18) for ωx , ωy into formula (20), we obtain the following
relation between the frequency ratio and the satellite mass ratio:
3q 2
N mD
= 2 − 4.
mC
p
In particular, positivity of the left-hand side entails
√
p
3
<
.
q
2
(22)
We remark in passing that this condition excludes the usual elliptic (or circular) oscillations corresponding to p = q = 1.
The frequencies ωx , ωy can then be expressed in terms of p and q:
pω0
ωx = p
,
q 2 − p2
qω0
ωy = p
.
q 2 − p2
10
4.1
Balanced formations avoiding collisions
We seek now formations of deputy satellites moving according to (21) subject to the following
conditions:
A. The arrangement of deputy satellites must satisfy the balance condition (19);
B. The satellites must never collide, i.e.
(xi (t), yi (t)) 6= (xj (t), yj (t))
for all t and i 6= j;
C. Optionally, we may wish to ensure that
(xi (t), yi (t)) 6= (0, 0)
for all i and t – in this case an additional satellite can be added at the center of the
“hub-and-spoke” system without collisions with this system.
We consider two types of uniform arrangement of deputy satellites. Type I is a uniform
arrangement of N satellites along a single Lissajous curve: the position of the i’th deputy
satellite is given by
i
i
xi (τ ) = x0 sin 2πp τ +
yi (τ ) = y0 sin 2πq τ +
+ ϕx ,
+ ϕy ,
N
N
where the phases ϕx , ϕy are the same for all satellites.
Type II is a uniform arrangement of N satellites along several Lissajous curves: the
position of the i’th deputy satellite is given by
i
i
xi (τ ) = x0 sin 2π pτ +
yi (τ ) = y0 sin 2π qτ +
+ ϕx ,
+ ϕy ,
N
N
where the phases ϕx , ϕy are the same for all satellites.
The following proposition summarizes properties of such formations with respect to the
above three conditions.
Proposition 1 Let N = 2, 3, . . . Denote ϕ0 =
qϕx −pϕy
.
π
a) For a Type I formation, the balance condition A is fulfilled if and only if neither p nor
q is divisible by N . For a Type II formation, the balance condition A is fulfilled for all
N.
b) For a Type I formation, the no-collision condition B is fulfilled iff the number ϕ0 +
(p − q)/2 is not an integer and N is co-prime with p and q.
11
c) For a Type II formation, in case N ≥ 3 the no-collision condition B is fulfilled if and
only if (ϕ0 + (p − q)/2) N is not divisible by the greatest common divisor of N and
q − p. In case N = 2 condition B is fulfilled iff ϕ0 is not an integer.
d) Lissajous curve (21) goes through the origin (0, 0) iff ϕ0 is an integer. It follows that
a Type I formation satisfies condition C iff ϕ0 is not an integer. A Type II formation
satisfies condition C iff ϕ0 is not of the form a + 2b(q − p)/N with integer a, b.
The proof of Proposition 1 is given in Appendix A. Note that, given p, q, admissible
formations (satisfying all three conditions A-C) exist for all N = 2, 3, . . . in case of Type II,
but not for all N in case of Type I.
In Table 1 we list parameters of admissible formations for all p, q ≤ 4 subject to condition
(22).
p/q(= ωx /ωy )
N mD /mC
1/2
1/3
2/3
1/4
3/4
8
23
11/4
44
4/3
Admissible N for Type I
3,
2,
5,
3,
5,
5,
4,
7,
5,
7,
7, . . .
5, . . .
11, . . .
7, . . .
11, . . .
Table 1: Admissible formations for small values of p, q.
In Fig. 2 we show several examples of formations of Types I and II satisfying conditions
A,B,C with p, q subject to condition (22).
4.2
Entanglement
Practical implementation of the introduced formations requires to resolve another issue. The
tethers’ ends are not attached to the main satellite at exactly the same point. Even if the
tethers are connected very close to each other, they still have nonzero thickness. Relative
motion of satellites in formations of Type I and Type II may cause not only contacts between
the tethers, but also tethers entanglement.
We will distinguish two kinds of entanglement: one that can be canceled out by rotating
the main satellite about the z axis as shown in Fig. 3a, and one that can not be eliminated by
such rotations (Fig. 3b). These two kinds will be referred to as weak and strong entanglement,
respectively.
We can define entanglement rigorously by applying simple topological concepts [2] to the
dynamics of the deputy satellites. Consider the motion of a deputy satellites formation as a
12
y
Type I, p=1, q=2, N=3, phi0=1/4
Type I, p=1, q=2, N=5, phi0=1/4
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−1.0
−1.0
−0.5
0.0
0.5
1.0
−1.0
Type I, p=1, q=3, N=2, phi0=1/2
y
Type II, p=1, q=2, N=2, phi0=1/2
1.0
−0.5
0.0
0.5
1.0
−1.0
Type I, p=3, q=4, N=5, phi0=1/4
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
−0.5
−0.5
−0.5
−1.0
−1.0
−0.5
0.0
0.5
1.0
x
0.0
0.5
1.0
Type II, p=2, q=3, N=2, phi0=1/2
1.0
−1.0
−0.5
−1.0
−1.0
−0.5
0.0
0.5
1.0
−1.0
−0.5
x
0.0
0.5
1.0
x
Figure 2: Examples of formations of Types I and II.
continuous map
M0 (τ ) =
x1 (τ )
y1 (τ )
xN (τ )
x2 (τ )
, τ ∈ [0, 1].
,...
,
yN (τ )
y2 (τ )
By periodicity, M0 (0) = M0 (1). Entanglement means the impossibility to “straighten” this
map subject to topological constraints. Precisely, let us say that M0 is strongly homotopic
to the constant map
xN (0)
x2 (0)
x1 (0)
, τ ∈ [0, 1] ,
,...
,
M1 (τ ) =
yN (0)
y2 (0)
y1 (0)
if there is a continuous family of maps Mα (τ ) , τ ∈ [0, 1] , α ∈ [0, 1] , that deforms M0 into
M1 without collisions and subject to the periodicity constraint
Mα (0) = Mα (1) , α ∈ [0, 1].
(23)
If M0 is not strongly homotopic to M1 , then we call M0 weakly entangled.
In order to define strong entanglement, we say that M0 is weakly homotopic to the
constant map M1 if the constraint (23) is replaced by the relaxed constraint
Mα (0) = Rϕ(α) Mα (1) , α ∈ [0, 1] ,
13
D2
D3
D2
D1
D2
D1
D3
C
D1
o
o
360 turn
of the D1D2D3 subsystem
D3
360 turn
of the main satellite C
C
D2
C
(a)
D2
D1
D3
D1
D3
C
C
(b)
Figure 3: Examples of weak and strong entanglement. In the case of weak entanglement
(a) tethers can be disentangled by rotation of the main satellite. Tethers in strong
entanglement (b) can not be disentangled by rotation of the main satellite.
where Rϕ(α) is the rotation of the satellite positions in the xy plane by continuously varying
angles ϕ(α) common to all deputy satellites, for some choice of the continuous function ϕ(α).
We then say that M0 is strongly entangled if M0 is not weakly homotopic to M1 .
The above definitions do not involve the main satellite. We assume that the initial state
is non-entangled, and the geometry of the tether attachment points to the main satellite is
reflected in the initial positions of the deputy satellites as shown in Fig. 3.
An obvious example of a weak, but not strong entanglement results from a circular motion
as shown in Fig. 3a. Also note that a motion of N = 2 deputy satellites is never strongly
entangled but may be weakly entangled.
By considering a few simple examples, it is easy to see that entanglement does take place
in some formations of Type I and Type II and does not in others. For instance, there is a
strong braid-like entanglement for the Type I formation with p = 1, q = 2, and N = 3 (top
left plot in Fig. 2), weak entanglement for the Type I formation with p = 1, q = 3, and
N = 2 (bottom left plot in Fig. 2), and no entanglement for the Type II formations with
N = 2 and p = 1, q = 2 or p = 2, q = 3, and N = 2 (top and bottom right plots in Fig. 2).
Complete analysis of entanglement in our formations appears to be relatively complex
mathematically, hence we will not attempt it here. We will, however, state a simple proposition involving pairwise relations between satellites. Namely, for any pair (i, j) of deputy
satellites consider the winding number wi,j defined as the number of turns that the j’th
satellite makes about the i’th satellite in the xy-plane over the period:
wi,j
1
=
arg
2π
xj (τ ) − xi (τ )
yj (τ ) − yi (τ )
1
τ =0
Clearly, wi,j = wj,i .The winding numbers are obviously invariant under a strong homotopy.
14
Under a weak homotopy, they are incremented by an amount common to all pairs and equal
1
to the number of full turns 2π
ϕ(α)|1α=0 . Since winding is absent in the straightened map
M1 , the numbers wi,j can serve to establish sufficient (but not necessary) conditions of
entanglement.
Proposition 2
a) Consider a formation of Type I or Type II moving without collisions. If at least one of
p, q is even, then wi,j = 0 for all pairs. Otherwise, wi,j = ±1 for all pairs.Consequently,
if both p, q are odd, then there is at least a weak entanglement.
b) Moreover, consider a formation of Type I and suppose that both p, q are odd so that
wi,j = ±1 by a). Then both values +1 and −1 are encountered among the winding
numbers wi,j iff neither of the numbers q−p, q+p is divisible by 2N . As a consequence,
if neither of q − p, q + p is divisible by 2N , then the entanglement is not only weak,
but also strong.
Proposition 2 is proved in Appendix B. It leads to an expected result if applied to a
circular or elliptic motion (p = q = 1): part a) confirms a weak entanglement, and part b)
ensures that all winding numbers are equal since p − q = 0, so the sufficient condition for a
strong entanglement is not fulfilled.
Note also that for odd p, q one of the numbers q − p, q + p is always divisible by 4, so
statement b) agrees with our earlier remark that there can be no strong entanglement for
N = 2.
4.3
Second order perturbation theory
In contrast to the general stability of the “hub-and-spoke” system established in Section
3.4 in the sense of smallness of deviation from the vertical equilibrium, we do not expect
the motion of deputy satellites along Lissajous curves to be stable. This motion is a subtle
phenomenon which, in particular, is strongly affected by nonlinearities and can be maintained
without additional control action only at relatively small oscillation amplitudes.
In this section we examine the nonlinearity effects by deriving second order corrections to
the evolution equations in case k/3ω02 mD ≫ 1 when one can neglect the tethers’ extensibility
and put l0 = l∗ . If we assume that the center of mass of the system rests at the origin, then
the configuration can be parametrized by the 2N coordinates xi , yi , i = 1, .., N, of the
15
auxiliary satellites. Specifically, the remaining coordinates are given by
N
N
mD X
mD X
xi , yC = −
yi ,
mC i=1
mC i=1
N q
mr X 2
zC =
l − (xi − xC )2 − (yi − yC )2 ≈
mC i=1
)
(
N
1 X
mr
2
2
,
N l0 −
(xi − xC ) + (yi − yC )
≈
mC
2l0 i=1
q
mr
zi = zC − l2 − (xi − xC )2 − (yi − yC )2 ≈ −
l0 +
mD
)
(
N
X
1
m
r
.
+
(xi − xC )2 + (yi − yC )2 −
(xk − xC )2 + (yk − yC )2
2l0
mC k=1
xC = −
Here and below we consistently keep terms only up to second order in xi , yi , ẋi , ẏi , ẍi , ÿi .
The x, y components Ti,x , Ti,y of the tension acting on the i’th deputy satellite can be
expressed through the z component by
Ti,x =
xi − xC
xi − xC
yi − yC
yi − yC
Ti,z ≈ −
Ti,z , Ti,y =
Ti,z ≈ −
Ti,z .
zi − zC
l0
zi − zC
l0
The z component Ti,z , in turn, is found from the corresponding HCW equation:
Ti,z
mr
≈ 2ω0 ẋi − 3ω02 zi ≈ 2ω0 ẋi + 3ω02
l0 .
m
mD
The remaining two HCW equations then yield
xi − xC
2 mr
ẍi ≈ 2ω0 żi −
2ω0 ẋi + 3ω0
l0 ,
l0
mD
yi − yC
2
2 mr
ÿi ≈ −ω0 yi −
2ω0 ẋi + 3ω0
l0 .
l0
mD
Differentiating zi and retaining our notation for ωx , ωy from (17)-(18), we finally obtain
2ω0
2
ẍi + ωx (xi − xC ) ≈
− ẋC (xi − xC ) + (ẏi − ẏC ) (yi − yC ) −
l0
N
mr X
(ẋj − ẋC ) (xj − xC ) + (ẏj − ẏC ) (yj − yC ) ,
−
mC j=1
ÿi + ωy2 yi − 3ω02
2ω0
mr
yC ≈ −
ẋi (yi − yC ) ,
mD
l0
where we have placed first order terms on the left and second order terms on the right.
16
The obtained equations can be used to find anharmonic corrections to a particular small
harmonic oscillation. The usual procedure is to substitute the harmonic oscillation in the
right-hand side and find the correction, to leading order, by solving the resulting nonhomogeneous linear equation [11]. If the right-hand side contains secular terms, i.e., those
whose frequencies match some of the eigenfrequencies of the linear equation, then, additionally, the solution’s frequencies need to be adjusted to prevent its non-physical growth.
In the case at hand we take the motion of a Type I or II formation with a small amplitude
as a base harmonic oscillation that we denote xi,0 (t) , yi,0 (t). From the balance condition
(19) we have xC (t) = yC (t) ≡ 0, so that the evolution equations simplify to
#
"
N
X
2ω
m
0
r
ẍi + ωx2 (xi − xC ) ≈
(ẋj,0 xj,0 + ẏj,0 yj,0 ) ,
ẏi,0 yi,0 −
l0
mC j=1
ÿi + ωy2 yi − 3ω02
2ω0
mr
yC ≈ −
ẋi,0 yi,0 .
mD
l0
The linear terms describe oscillations with frequencies
√
ωx , ωy , ωCx = 3ω0 , ωCy = 2ω0 ,
where the latter two correspond to the motion of the center of mass of auxiliary satellites or,
equivalently, to the motion of the main satellite. The terms ẋk,0 xk,0 , ẏk,0 yk,0 , and ẋi,0 yi,0 on
the right-hand side result in oscillations with frequencies 2ωx , 2ωy , and ωy ± ωx , respectively.
Though the second order correction of motion is obviously present for each deputy satellite
for any choice of system parameters, it is possible to choose parameters so as to make the
second order correction completely vanish for the main satellite. The evolution equations for
the main satellite are derived by adding up the equations for the deputy satellites:
!
N
N
X
X
2ω0
N mD
mC
mC
2
ẏj,0 yj,0 −
ẋj,0 xj,0 ,
−
ẍC + ωCx xC ≈
mD
l0
mC + N mD j=1
mC + N mD j=1
−
N
2ω0 X
mC
2
ÿC + ωCy
yC ≈ −
ẋj,0 yj,0 .
mD
l0 j=1
Note that the N mD tends to be much larger than mC for a system satisfying our assumptions (see Table 1), so the motion of the main satellite tends to be generally affected by
second order corrections much stronger than the motion of the deputy satellites. However,
contributions from different j’s here may cancel out.
Proposition 3 Consider a Type I formation such that neither of 2p, 2q, q ± p is divisible
by N . Then
N
X
k=1
ẋk,0 (t)xk,0 (t) =
N
X
ẏk,0 (t)yk,0 (t) =
k=1
N
X
k=1
ẋk,0 (t)yk,0 (t) ≡ 0
for all t, so that the right-hand sides of the above equations vanish identically.
17
The proof is elementary, and we omit it. Examples of parameter sets fulfilling the proposition’s hypothesis are N = 5, p = 1, q = 2 and N = 5, p = 3, q = 4. Our numerical
experiments below confirm that in these cases the main satellite is indeed stable in contrast
to the generic settings at the same oscillation amplitude of deputy satellites.
5
5.1
Examples and numerical simulations
Setup and general observations
In our numerical simulations we consider oscillations of the system with rigid tether at small
angles and with equal amplitudes in x and y:
x0 = y0 = a = κrad l0
Here a is the linear amplitude and κrad is the corresponding angular amplitude expressed
in radians; the same angular amplitude expressed in degrees is denoted κdeg . In all our
experiments the dimensionless coefficient k/3ω02 mD characterizing tether rigidity falls in
the range [3 · 102 , 103 ] and κdeg ≤ 6◦ (to justify the above linear relation between a and
κrad ). In order to numerically examine the stability of the Lissajous motion we introduce
quantities characterizing relative deviations of the main and deputy satellites from their
theoretical positions obtained in the linear approximation. Specifically, we consider the
relative deviation of the main satellite’s numerically computed position from the z axis
1q 2
2
δC (t) =
(t),
xC,num (t) + yC,num
a
and the mean relative deviation of the deputy satellites’ numerically found trajectories from
the theoretical Lissajous curves:
N
1 X
δD (t) =
N a i=1
q
(xi,num (t) − xi,Liss (t))2 + (yi,num (t) − yi,Liss (t))2
These quantities can be compared with the minimum distance between different deputy
satellites on their theoretical trajectories:
q
1
δmin = min (xi,Liss (t) − xj,Liss (t))2 + (yi,Liss (t) − yj,Liss (t))2
a t,i6=j
Here we neglect the difference in z coordinates of the deputy satellites since for small amplitude oscillations this difference has higher order of smallness.
If each satellite remains at all times within a distance of δmin /2 from its theoretical
position on the Lissajous curve, then all satellites are guaranteed to avoide collisions. We
will consider a slightly relaxed condition of stability
δD (t) <
18
δmin
2
(24)
that constrains only the mean deviation δD (t) of the deputy satellites. The time interval
during which this condition holds can be roughly considered as a “system stability interval”.
We restrict ourselves to the simplest ratio ωx : ωy = 1 : 2 achieved at p = 1, q = 2. In
this case, the theoretical mass ratio is
N mD
= 8,
mC
and the full period TL of system oscillations is related to the orbital period T0 = 1/(2πω0 )
by
p
√
TL = q 2 − p2 T0 = 3T0 ,
so that the non-dimensional time
t
.
3T0
We consider the three different configurations of the system shown in the first row of
Fig. 2. For each configuration, we perform simulations for κdeg = 1◦ and κdeg = 3◦ for 10
orbital periods. In our simulations, the central satellite moves along a geostationary orbit,
and tether length l0 = 10000 m.
δ
τ=√
0.30
0.25
0.20
0.15
0.10
0.05
0.00
0
0.5
δ
0.4
ϰdeg =1 ◦
δC
δD
δmin/2
2
4
6
8
ϰdeg =3 ◦
δC
δD
δmin/2
0.3
0.2
0.1
0.0
0
10
2
4
orbit number
6
8
10
Figure 4: Numerically observed relative deviations of the satellites from the theoretical
trajectories for the Type I formation with N = 3 deputy satellites. In this case the
minimum relative distance between deputy satellites on the theoretical trajectories is
δmin = 0.60. At κdeg = 1◦ stability condition (24) holds with a large margin for the
whole simulation interval of 10 orbital periods, while at κdeg = 3◦ it breaks down after
six orbital periods.
19
The results for Type I formation at N = 3 (Fig. 4) show that the deviations of the
main satellite are initially much larger than those of the deputy satellites. This is not
surprising, since the main satellite is eight times lighter. However, deviations of the deputy
satellites approximately linearly accumulate with time, and eventually catch up with those
of the main satellite. In case κdeg = 3◦ relative deviations are much larger than in case
κdeg = 1◦ : approximately three times larger for δC and six times larger for δD . The second
order perturbation theory in Section 4.3 suggests a linear dependence of relative deviations
on κdeg , but for κdeg = 3◦ deviations of the main satellite from the equilibrium position
are already comparable to the oscillation amplitude, so this perturbation theory is not truly
applicable here. In case κdeg = 1◦ the stability condition (24) holds with a large margin for
the whole simulation interval, while at κdeg = 3◦ it breaks down after six orbital periods.
0.20
δ
0.15
ϰdeg =1 ◦
δC
δD
δmin/2
0.10
0.05
0.00
0
0.20
δ
0.15
2
4
6
8
ϰdeg =3 ◦
δC
δD
δmin/2
0.10
0.05
0.00
0
10
2
4
orbit number
6
8
10
Figure 5: Deviations for Type I and N = 5. In this case δmin = 0.43. Stability condition
(24) holds with a large margin for the whole simulation interval of 10 orbital periods
for both κdeg = 1◦ and κdeg = 3◦ .
The results for Type I formation at N = 5 (Fig. 5) are drastically different due to
the cancellation of second order corrections of the main satellite’s motion pointed out in
Proposition 3. Not only is the main satellite almost immobile, but also the deviations of the
deputy satellites are from two to four times smaller than in the previous case. In particular,
the no-collision condition (24) holds at κdeg = 3◦ throughout the whole simulation interval.
Note also that in the κdeg = 3◦ case the deviations of deputy satellites are approximately
three times as large as in the κdeg = 1◦ case, in good agreement with the second order
perturbation theory which is now applicable since the deviations are small.
Results in the case of Type II formation at N = 2 (Fig. 6) are on the whole similar to
20
δ
0.6
0.5
0.4
0.3
0.2
0.1
0.0
0
1.0
δ
0.8
ϰdeg =1 ◦
δC
δD
δmin/2
2
4
6
8
ϰdeg =3 ◦
0.6
δC
δD
δmin/2
0.4
0.2
0.0
0
10
2
4
orbit number
6
8
10
Figure 6: Deviations for Type II and N = 2. In this case δmin = 1.32. At κdeg = 1◦
stability condition (24) holds with a large margin for the whole simulation interval of
10 orbital periods, while at κdeg = 3◦ it breaks down after six orbital periods.
those obtained in the first case (Type I, N = 3). The deviations are in fact now higher than
in that case, but this is somewhat compensated by the larger δmin , so that the no-collision
condition at κdeg = 3◦ again holds about up to the half of the simulation interval.
5.2
Frequency adjustment by tuning the mass ratio
The specific pattern of deviation growth observed in the above examples (linear growth in
time with superimposed periodicity) strongly suggests that this growth is largely due to the
gradual shift of the trajectories occurring because the ratio of the true x and y system’s
frequencies does not exactly match the approximate value α0 = N mD /mC = 8. We can
expect to negate these shifts by adjusting the mass ratio:
N mD
= α0 − ∆α.
mC
We find the appropriate ∆α numerically, by minimizing the maximum deviation of deputy
satellites for 10 orbital periods:
max δD (t) → min
∆α
t∈[0,10T0 ]
We perform these optimizations at different angles κdeg for the second, most stable configuration from the previous section (Type I, N = 5). The results are shown in Table 2. The
21
κdeg
1◦
2◦
3◦
4◦
5◦
6◦
∆αopt
0.021
0.039
0.071
0.116
0.169
0.230
max δD (t)
t∈[0,10T0 ]
max δD (t)
t∈[0,10T0 ]
max δD (t)
t∈[0,30T0 ]
w/o adjustment
with adjustment
with adjustment
0.0408
0.0766
0.1230
0.1780
0.2600
0.3700
0.0210
0.0368
0.0541
0.0702
0.0874
0.1110
0.0375
0.0530
0.0814
0.1280
0.1650
0.2020
Table 2: Results of the mass ratio adjustment experiments. At each angle κdeg we
numerically find the optimal adjustment ∆αopt . The maximum deviations of the main
satellite from the vertical without adjustment or with the optimized adjustment are
shown in the third and fourth columns, respectively. The last column shows results of
simulations with adjustment spanning 30 orbital periods.
results clearly show a big improvement over the earlier results obtained without mass ratio
adjustment. In Fig. 7 the deviations of the adjusted systems with different values of κdeg
are plotted for 30 orbital periods.
6
Conclusions
We found that subject to appropriate choice of “hub-and-spoke” system parameters deputy
satellites can move along Lissajous curves so that the system stays in free motion (i.e. no
fuel is consumed in the nominal operation mode). The main satellite is in a state of relative
equilibrium on the local vertical, passing through the system CoM.
Our analysis shows the existence of rather nontrivial relations between the system’s
parameters that, when satisfied, produce a well-balanced system without collisions between
the deputy satellites or tethers. Certain configurations determined by the proper choice of
parameters allow placing an additional satellite at the center of the “hub-and-spoke” system
without collisions with the other parts of the system.
Of course, the intricate way in which the deputy satellites change their positions is challenging for technical implementation. In particular, the tethers must be attached to the main
satellite so as not to intertwine. Analysis of the tether entanglement in terms of homotopy
provides another set of constraints that the system parameters should satisfy.
One of the most curious results in this paper is the study of nonlinear effects, which
allows to formulate yet another set of conditions for the systems parameters to cancel out
the second order corrections in the systems’ equations of motion. Our numerical experiments
with two parameter sets (for 5 deputy satellites) satisfying these conditions corroborate the
theoretically predicted system’s stability
The application of the proposed motion pattern is limited to small angular deviations of
22
0.20
δ
0.15
ϰdeg =1 ◦
δC
δD
δmin/2
0.10
0.05
0.00
0
0.20
δ
0.15
5
10
15
20
25
ϰdeg =3 ◦
δC
δD
δmin/2
0.10
0.05
0.00
0
0.20
δ
0.15
5
10
15
20
25
30
ϰdeg =6 ◦
δC
δD
δmin/2
0.10
0.05
0.00
0
30
5
10
15
orbit number
20
25
30
Figure 7: Deviations in the Type I, N = 5 formation with the optimally adjusted
satellite mass ratio. Stability condition (24) holds for the whole simulation interval of
30 orbital periods for all three considered values of κdeg .
the tethers from the local vertical. The possibility to extend it to large deviations from the
vertical could be a subject of further investigation.
Acknowledgment
The authors first conceived the idea of the motion described in this paper during the dynamical analysis of the rotating multi-tethered satellite system [1], and the authors sincerely
thank their collaborators Didier Alary, Kirill Andreev, Pavel Boyko, Elena Ivanova, and
Cyrille Tourneur for the warm and stimulating atmosphere of that study. The work of one
of the authors (DY) on the present paper was supported by Russian Science Foundation
(project 14-50-00150).
23
Appendix A
Proof of Proposition 1
Statements a) and d) are very simple, so we only provide proofs for b) and c).
b) Suppose that the satellite i collides with the satellite j at the time moment τ . The
condition xi (τ ) = xj (τ ) admits two series of solutions:
1. 2πp τ + Ni + ϕx = 2π p τ + Nj + nx + ϕx with some integer nx , that is
i−j
nx
=
.
N
p
2. 2πp τ +
i
N
+ ϕx = 2π
+ nx − p τ + Nj − ϕx with some integer nx , that is
i+j
ϕx
1 1 1
−
.
+ nx −
τ=
2 p 2
π
N
1
2
The condition yi (τ ) = yj (τ ) has similar series but with q instead of p, ϕy instead of ϕx
and ny instead of nx .
Thus there are four possibilities for the satellites to collide: an element in series 1 or
2 for the coordinate x must occur simultaneously with an element in series 1 or 2 for the
coordinate y. Let us deal with these cases one by one.
1. Series 1 for x and series 1 for y – impossible because p and q are co-prime and −N <
i − j < N.
2. Series 2 for x and series 1 for y. Series 1 for y has non-trivial solutions iff N and q are
not co-prime. If a non-trivial solution exists, the time moment τ is found from series
2 for x. Thus in this series collisions occur iff N and q are not co-prime.
3. Series 1 for x and series 2 for y. Similarly, collisions occur iff N and p are not co-prime.
4. Series 2 for x and series 2 for y. Equating τ from both series, we obtain
1 1
1 1
ϕx
ϕy
=
+ nx −
+ ny −
p 2
π
q 2
π
whence
nx q−ny p=
qϕx − pϕy p − q
+
.
π
2
Since p and q are co-prime, with nx , ny running over all possible integers, the left side
also runs over all possible integers, i.e. collisions in this case happen iff the right side is an
integer.
c) Like in the proof of b), we obtain two series of relations from the condition xi (τ ) =xj (τ ).
However, the first series takes the form
i−j
=nx
N
24
and has no non-trivial solutions, because −N <i − j<N . This leaves the second series, which
has the form
1 1
ϕx i + j
τ=
.
+ nx −
−
2p 2
π
N
Like before, we equate τ from the series for x and y and obtain
1 1
ϕx i + j
ϕy i + j
1 1
=
+ nx −
−
+ ny −
−
p 2
π
N
q 2
π
N
whence
N (qnx − pny ) − (q − p) (i + j) =
qϕx − pϕy p − q
+
π
2
N.
As nx , ny run over all integer values, the expression qnx −pny also takes all integer values.
If N ≥ 3, the sum i + j for all possible pairs of different numbers from 1 to N assumes all
possible integer values modulo N . Thus the expression in the left side takes all possible
integer values divisible by the greatest common divisor of N and q − p. On the other hand,
if N = 2, then i + j = 3 and the above relation simplifies to the requirement that
qϕx − pϕy
π
be integer.
Appendix B
Proof of Proposition 2
a) For a closed curve (x (τ ) , y(τ ))τ ∈[0,1] not containing the origin, the number w of its rotations about the origin can be written as
w=
1 X
sign (x (τ ) ẏ (τ )),
2
τ :y(τ )=0
assuming that the intersections of the curve with the x coordinate axis are non-degenerate.
We will apply this formula to
xij (τ ) = xj (τ ) − xi (τ ) , yij (τ ) = yj (τ ) − yi (τ ) .
We consider separately the two types of formations.
Type I. In this case
πp(j − i)
j+i
xij (τ ) = 2x0 sin
+ ϕx ,
cos 2πp τ +
N
2N
j+i
πq(j − i)
+ ϕy ,
cos 2πq τ +
yij (τ ) = 2y0 sin
N
2N
25
so that
wi,j
πq(j − i)
πp(j − i)
sign sin
w∗ = ±w∗ ,
= sign (x0 y0 ) sign sin
N
N
where w∗ is the winding number of the Lissajous curve
x∗ (τ ) = cos (2πpτ + ϕx ) , y∗ (τ ) = cos (2πqτ + ϕy ) .
We have
1
y∗ (τ ) = 0 for τ =
2q
ϕy 1
s−
−
π
2
, s = 1, ..2q.
Moreover, with this choice of τ we have ẏ∗ (τ ) > 0 for even s and ẏ∗ (τ ) < 0 for odd s.
Applying the formula for the winding number, we obtain
2q
1X
πps
s
+ ϕ∗
(−1) sign cos
w∗ =
2 s=1
q
with some (unimportant) phase constant ϕ∗ .
Now consider separately the cases when both p, q are odd and when one of them is even.
Let both p, q be odd. Since p, q are co-prime, the values (πps/q)mod 2π run
over the
2q−1
πps
values in the set Z2q = {πr/q}r=0 as s runs over 1, 2, .., 2q. Then the quantity cos q + ϕ∗
takes equally many positive and negative values as s runs over 1, 2, .., 2q. Now, if s runs only
over even values 2, .., 2q, then (πps/q)mod 2π runs over the values in the set Zq ={2πr/q}q−1
r=0
.
Then, since q is odd, the number of positive and negative values taken by cos
as s runs only even values 2, .., 2q, differs by 1. It follows that
w∗ =
πps
q
+ ϕ∗ ,
1
(±1 − (∓1)) = ±1.
2
Now let one of p, q be even; without loss of generality we can assume that it is
q. We can
πps
then repeat the above argument, but since q is now even, we conclude that cos q + ϕ∗
takes equally many negative and positive values as s runs over the even or odd subset of
1, 2.., 2q. It follows that w∗ = 21 (0 − 0) = 0.
We have thus proved statement a) for Type I formation.
Type II. In this case
π(j − i)
j+i
xij (τ ) = 2x0 sin
+ ϕx ,
cos 2π pτ +
N
2N
j+i
π(j − i)
cos 2π qτ +
+ ϕy ,
yij (τ ) = 2y0 sin
N
2N
which leads to
wi,j
2q
sign (x0 y0 ) X
ps i + j q − p
s
(−1) sign cos π
=
+
+ ϕ∗
2
q
N
q
s=1
26
with some phase constant ϕ∗ not depending on i, j. The statement a) then follows just like
in case of Type I.
b) From the proof of a), wi,j = w∗ sign (x0 y0 ) sign(fj−i ), where we set
πps
πqs
fs = sin
sin
.
N
N
Note first that if q − p is divisible by 2N , then fs = sin2 (πps/N ) > 0 for all s not divisible
by N (recall that N and p are co-prime by the assumption of no collisions). Accordingly,
wi,j = w∗ sign (x0 y0 ) for all i, j = 1, .., N with i 6= j. Similarly, wi,j = −w∗ sign (x0 y0 ) if q + p
is divisible by 2N .
Now suppose that neither of q − p, q + p is divisible by 2N . We need to show that fs
takes both positive and negative values
as s runs over 1, 2, .., N − 1. Since fN = f0 = 0 and
P2N
f2N −s = fs , it suffices to show that s=0−1 fs = 0. But that immediately follows from the
hypothesis and the identity
πs(p − q)
πs(p + q)
1
cos
.
− cos
fs =
2
N
N
The presence of strong entanglement follows from the presence of different winding numbers,
since under a weak homotopy the winding number (possibly nonzero) must be the same for
all pairs.
References
[1] Alary, D., Andreev, K., Boyko, P., Ivanova, E., Pritykin, D., Sidorenko, V., Tourneur,
C., Yarotsky, D.: Dynamics of multi-tethered pyramydal satellite formation. Acta Astronautica. 117, 222-232 (2015)
[2] Armstrong, M.A.: Basic topology. Springer, New York (1983)
[3] Arnold, V.I.: Ordinary differential equations. MIT Press, Cambridge, Massachussets
(1978)
[4] Avanzini, G., Fedi, M.: Refined dynamical analysis of multi-tethered satellite formations. Acta Astronautica. 84, 36-48 (2013).
[5] Avanzini, G., Fedi, M.: Effects of eccentricity of the reference orbit on multi-tethered
satellite formations. Acta Astronautica. 94, 338-350 (2014).
[6] Bekey, I.: Tethers open new space options. Astronautics and Aeronautics, 21, 23-40
(1983).
[7] Cai, Z., Zhao, J., Peng, H., Qi, Z.: Nonlinear control of rotating multi-tethered formations in halo orbits. Int. J. Comput. Methods, 11, 1344008 [26 pages](2014).
27
[8] Celletti, A., Sidorenko, V.V.: Some properties of dumbbell satellite attitude dynamics.
Cel. Mech. Dyn. Astron., 101, 105-126 (2008).
[9] Clohessy, W.H., Wiltshire, R.S.: Terminal guidance system for satellite rendezvous. J.
Aerosp. Sci., 27, 653-658 (1960).
[10] Gantmacher, F.: Applications of the theory of matrices, Interscience, New York (1959)
[11] Nayfeh, A.H.: Perturbations methods, Wiley-Interscience, New York (1973)
[12] Panosian, S., Seubert, C.R., Schaub, H.: Tethered Coulomb Structure Applied to Close
Proximity Situational Awareness AIAA Journal of Spacecraft and Rockets, 49, 11831193 (2012).
[13] Pizarro-Chong, A., Misra, A.K.: Dynamics of multi-tethered satellite formations containing a parent body. Acta Astronautica, 63, 1188-1202 (2008).
[14] Seubert, C.R., Schaub, H.: Tethered Coulomb Structures: Prospects and Challenges,
Journal of Astronautical Sciences, 57, 347-348 (2009)
[15] Seubert, C.R., Panosian, S., Schaub, H.: Attitude and Power Analysis of MultiTethered, Two-Node Tethered Coulomb Structures AIAA Journal of Spacecraft and
Rockets, 48, 1033-1045 (2011)
[16] Sidorenko, V.V., Celletti, A.: A “spring-mass” model of tethered satellite systems:
properties of planar periodic motions. Cel. Mech. Dyn. Astron., 107, 209-231 (2010).
[17] Wong, B., Misra, A.: Planar dynamics of variable length multi-tethered spacecraft near
collinear Lagrangian points. Acta Astronautica, 63, 1178-1187 (2008).
[18] Zhao, J., Cai, Z.: Nonlinear dynamics and simulation of multi-tethered satellite formations in Halo orbits. Acta Astronautica, 63, 673-681 (2008).
28