Analysis on reachable set for spacecraft relative motion under low-thrust

Automatica 115 (2020) 108864

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Brief paper

Analysis on reachable set for spacecraft relative motion under low-thrust✩ Hao Sun a , Junfeng Li b , a b

∗

Yau Mathematical Sciences Center, Tsinghua University, 100084, Beijing, People’s Republic of China Lab of Dynamics and Control, Tsinghua University, 100084, Beijing, People’s Republic of China

article

info

Article history: Received 23 March 2019 Received in revised form 10 November 2019 Accepted 11 January 2020 Available online xxxx Keywords: Reachable set Hill’s equations Spacecraft relative motion Energy-limited constraint

a b s t r a c t Hill’s equations describe the linearized relative motion of deputy spacecraft with respect to the chief spacecraft which moves around the Earth in a circular orbit. The reachable set is defined as the set of all possible final positions of the deputy spacecraft that can be reached by a control force under certain constraints. The computation of the reachable set has been a focus of recent study, and different numerical methods are proposed. In this article, we use rigorous mathematical methods to study the structure and properties of the reachable set under the energy-limited constraint. We give an analytical expression for the computation of reachable set by using the methods of Hilbert space theory. According to the obtained analytical expression, the shape, volume, surface area and other properties of the reachable set are studied. In the end, the conclusions are verified numerically. © 2020 Published by Elsevier Ltd.

1. Introduction Since space missions are becoming increasingly complex, more than one spacecraft is used to complete a mission cooperatively. It is often required that one spacecraft has a given relative position with respect to another spacecraft, and the control of spacecraft relative motion has been a focus of research work. The use of more than one spacecraft in a flight mission also introduces fundamental research problems such as formation keeping, spacecraft docking, and collision avoidance (Fehse, 2003; Schlanbusch, Kristiansen, & Nicklasson, 2011). All the problems are closely related to the fundamental study on spacecraft relative motion. Hill’s equations describe the linearized relative motion of deputy spacecraft with respect to the chief spacecraft which moves around the Earth in a circular orbit. Generally, the initial state of deputy spacecraft is known and the control force is required to satisfy the energy-limited constraint. The reachable set is defined as the set of all possible final positions of the deputy spacecraft that can be reached by a control force under the constraint of energy. The reachable set reflects the operational capabilities and limitations of the spacecraft, which is crucial in ✩ This work was partially supported by NSFC, People’s Republic of China grants 11871299 and 11672146. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Denis Arzelier under the direction of Editor Richard Middleton. ∗ Corresponding author. E-mail addresses: [email protected] (H. Sun), [email protected] (J. Li). https://doi.org/10.1016/j.automatica.2020.108864 0005-1098/© 2020 Published by Elsevier Ltd.

flight missions as it can be used for (i) mission planning and evaluation of the feasibility of missions, and (ii) collision probability monitoring with the objective of collision avoidance. From this perspective, plenty of work has been performed on the study of the reachable set. In control literature, the optimal control theory can be applied to study the reachable set, which requires finding the solution of Hamilton–Jacobi equation (Blanchini, 1999). Kurzhanski and Varaiya (2000) developed the ellipsoidal approximation method for the computation of the reachable set. Chen, Erika, and Sriram (2013) described the tool Flow∗ and applied it to study the behavior of hybrid systems. Frehse, Guernic, et al. (2011) presented a scalable reachability algorithm for hybrid systems with non-deterministic dynamics. Xue, Li, and Baoyin (2015) proposed a method which computes the approximated upper bound on the reachable set of a spacecraft with a single impulse. Wen and Gurfil (2016) developed a method of computing the relative reachable set for a spacecraft when initial state uncertainties are given. Lee and Hwang (2018) proposed a systematical way to compute both the inner and outer ellipsoids which approximate the reachable set under the constraint of energy. In respect of the computation for the reachable set, only numerical methods have been proposed, which have significant computational loads. So far, there has been no pure analysis and the structure of the reachable set is unknown. In this article, we use rigorous mathematical methods to study the structure of the reachable set under the energy-limited constraint. The methods of Hilbert space theory are applied to give an analytical expression for computing the reachable set. According to the obtained analytical expression, we study the shape, volume, surface

2

H. Sun and J. Li / Automatica 115 (2020) 108864

Suppose the initial state of the deputy spacecraft is r (0) = r0 and )⊤ ( )⊤ ( and F = 0, 0, 0, f ⊤ , then r˙ (0) = v0 . Denote P = r ⊤ , (˙r )⊤ the state–space equation associated to Eq. (2) is

{

P˙ = AP + F P (0) = P0

,

(4)

where

( ) r0

P0 = Fig. 1. The relative motion of the deputy spacecraft with respect to the chief spacecraft.

area and other properties of the reachable set. In the end, we use numerical methods to verify the theoretical conclusions. The numerical results coincide well with our analysis. In many space missions, the low-thrust propulsion system has been widely used due to its high specific impulse and high fuel efficiency. The energy-limited constraint is frequently considered for a spacecraft with the low-thrust propulsion system. The computation of reachable set has been a focus of research work, as the reachable set reflects the operational capabilities and limitations of the spacecraft. The reachable set can also be used for monitoring collision probability and evaluating mission feasibility. This article thoroughly solves the computation problem of the reachable set under the energy-limited constraint. The explicit expressions obtained in this article can be used to calculate the accurate reachable set without any computational loads. This enables the spacecraft to maintain high-level and immediate cognitive awareness about its limitation and environment (neighboring spacecrafts). 2. Hill’s equations

r¨ = A1 r + A2 r˙ + f ,

⎞

0 0 ⎠ , A2 =

−ω2

)

.

(5)

P (t ) = e

At

(

t

∫

P0 +

−Aτ

e

F (τ ) dτ

)

.

(6)

0

3. Formulation of the control problem In this section, we give the formulation of Hill’s equations as a control problem. Suppose the initial state of the deputy spacecraft is (P0 ∈ R6 . The force is f ∈ ) control time is T (and the control ) L2 [0, T ] , R3 . The Hilbert space L2 [0, T ] , R3 is equipped with the inner product

⟨f , g ⟩L2 =

T

∫

f · gdτ ,

(7)

0

where · is the inner product of vectors in R3 . Denote U as the set of all admissible control forces and Pf (t ) is the solution (6) under the control force f ∈ U. The reachable set is defined as

{

} (I3 , 0) Pf (T ) |f ∈ U .

(8)

In this article, we consider the constraint of energy. The energy of the control force is 1

T

∫

2

|f (τ )|2 dτ = 0

1

∥f ∥2L2 ,

2

(9)

3 where | · | is the norm √is under √ of a vector√in R . The control force the constraint (i) 2E ≤ R, (ii) 2E = R, or (iii) r ≤ 2E ≤ R, where R, r > 0 are given constants. The set of admissible control forces is

U1 = f ∈ L2 [0, T ] , R3 |∥f ∥L2 ≤ R

{

(

)

}

. } | r ≤ ∥f ∥L2 ≤ R

U2 = f ∈ L2 [0, T ] , R3 |∥f ∥L2 = R

{

(

{

( 2

U3 = f ∈ L

)

[0, T ] , R

) 3

}

(10)

We define the map G:L

2

[0, T ] , R

(

3

)

→ R , f ↦→ (I3 , 0) e 3

AT

T

∫

e−Aτ F (τ ) dτ ,

(11)

0

then the reachable set is (I3 , 0) eAT P0 + G (Ui ) respectively, i = 1, 2, 3. By the linearity of the control system (Eq. (4)), we only need to consider the condition P0 = 0. In the later sections, we suppose the initial state of the deputy spacecraft is P0 = 0.

In this section, we use the ) of Hilbert Space Theory ( methods to study the map G : L2 [0, T ] , R3 → R3 , which leads to the structure and properties of the reachable set. We have the following theorem.

(2)

Theorem 1. The map G defined by the expression (11) is

where 0 0 0

I3 A2

(1)

Hill’s equations can be expressed in a more concise form:

3ω2 ⎝ 0 A1 = 0

0 A1

4. Application of Hilbert space theory

⎧ 2 ⎪ ⎨ x¨ − 2ω˙y − 3ω x = f1 y¨ + 2ω˙x = f2 . ⎪ ⎩ 2 z¨ + ω z = f3

⎛

(

By the Duhamel formula in differential equation theory, the solution of Eq. (4) is

E=

A brief introduction to Hill’s equations is given in this section. Hill’s equations describe the linearized relative motion of the deputy spacecraft with respect to the chief spacecraft, which orbits a planet with the central spherical gravity field (the Earth) in a circular orbit (see Fig. 1). The chief and deputy spacecrafts are considered as particles. Suppose O1 XYZ is an inertial frame and O1 is the center of the Earth. The chief spacecraft moves in a circular orbit centered at O1 in the XY -plane, and the angular velocity ω > 0 is a constant. The moving frame O2 xyz is defined as follows. At any moment, the point O2 is the position of chief spacecraft. The directions of O2 x, O2 z are the same as those of O1 O2 , O1 Z respectively. The direction of O2 y is the same as that of the velocity of chief spacecraft (perpendicular to O1 O2 ). The relative position of the deputy spacecraft in O2 xyz is r (t ) = (x (t ) , y (t ) , z (t ))⊤ , where ⊤ represents the transpose of a vector or a matrix. The control force is denoted as f (t ) = (f1 (t ) , f2 (t ) , f3 (t ))⊤ , and Hill’s equations are (Hill, 1878)

v0

,A =

(

0 −2ω 0

2ω 0 0

0 0 . 0

)

(3)

(i) (ii) (iii) (iv)

linear, bounded, surjective, Lipschitz continuous.

H. Sun and J. Li / Automatica 115 (2020) 108864

Proof. ( ) (i) For ∀λ1 , λ2 ∈ R and ∀f1 , f2 ∈ L2 [0, T ] , R3 , we have G (λ1 f1 + λ2 f2 ) = λ1 G (f1 ) + λ2 G (f2 ). Thus G is linear. (ii) Denote ∥ · ∥2 as the 2-norm of a matrix ( and | · | )as the Euclidean norm of a vector. For any f ∈ L2 [0, T ] , R3 , we have

T

J1 J= 0 0

|F (τ )| dτ ≤

T

dτ · 0

0

|f (τ )|2 dτ

) 12

0 ⎜1 ) ⎜ 0 ⎜0 0 , Q =⎜ ⎜0 J3 ⎝0 0

−2/(3ω)

−i/2

1 0 0 1 0

1 0

ω/2 iω 0

i/2 1 0 ω/2 −iω 0

0 0 1 0 0 iω

⎞

0 0 ⎟ ⎟ 1 ⎟ ⎟, 0 ⎟ 0 ⎠ −iω (16)

( J1 =

0 0

)

1 , J2 = 0

iω 0

(

)

0 , J3 = iω

( −iω

)

0 . −iω

0

(17)

Further, we have eAt = QeJt Q −1 and

⎛

T

∫

0 J2 0

(12)

0

(∫

In this section, we use the Gram–Schmidt orthogonalization to construct the orthonormal basis of (Ker (G))⊥ . In Eq. (6), the matrix A has eigenvalues {0, ±iω}. To calculate eAt , we can utilize the Jordan canonical form A = QJQ −1 , where

(

eJ1 t Jt e =⎝ 0 0

By the Hölder inequality, we have

∫

5. Construction of the orthonormal basis

⎛

|G (f )| ⏐∫ ⏐   ⏐ T −Aτ ⏐ ≤ eAT 2 ⏐⏐ e F (τ ) dτ ⏐⏐ 0 ∫ T  AT  ⏐ −Aτ ⏐ ⏐e F (τ )⏐ dτ ≤ e 2 . 0 ∫ T  AT   −Aτ  e  · |F (τ )| dτ ≤ e 2 2 0 ∫ T  AT   −Aτ  |F (τ )| dτ ≤ e 2 · sup e 2 · τ ∈[0,T ]

3

.

(13) eJ1 t =

0

Combining the inequalities (12) and (13), we have

    √ |G (f )| ≤ ∥f ∥L2 · eAT 2 · T sup e−Aτ 2 , τ ∈[0,T ]

(14)

thus G is bounded. (iii) The proof is constructive. We need to prove )⊤ deputy ( that the ∈ R3 . spacecraft can reach any final position r = x , y , z f f f f ( Denote) f1 (t ) = 2rf /T 2 and f2 (t ) = −4yf ωt − 3xf ω2 t 2 , 4xf ωt , zf ω2 t 2 /T 2 . It can be verified that r (t ) = t 2 rf /T 2 is the solution of Hill’s equations when the control force is f = f1 + f2 . Clearly, G (f ) = r (T ) = rf and map. ( G is a surjective ) (iv) For ∀f1 , f2 ∈ L2 [0, T ] , R3 , we have

  √   |G (f1 )− G (f2 )| ≤ ∥f1 − f2 ∥L2 eAT 2 · T sup e−Aτ 2 , τ ∈[0,T ]

(15)

thus G is Lipschitz continuous. Due to the continuity of( G, the set) Ker (G) = G−1 ({0}) is a closed linear subspace of L2 [0, T ] , R3 . Denote (Ker (G))⊥ as the orthogonal complement space of Ker (G), which is also a closed linear subspace. By( the orthogonal projection theorem (Yosida, ) 1978), we have L2 [0, T ] , R3 = Ker (G) ⊕ (Ker (G))⊥ , where ⊕ is the direct sum of subspaces. We restrict G to (Ker (G))⊥ and obtain the map GE : (Ker (G))⊥ → R3 , which is also linear, bounded, and continuous. We have the following theorem. Theorem 2. (i) The map GE is an isomorphism and also a homeomorphism. (ii) (Ker (G))⊥ is a three-dimensional linear space. Proof. ( ) (i) For ∀x ∈ L2 [0, T ] , R3 , suppose x = x1 + x2 and x1 ∈ Ker (G) , x2 ∈ (Ker (G))⊥ . Then we have G (x) = G (x1 ) + G (x2 ) = G (x2 ). Since G is surjective, GE is also surjective. Suppose x1 , x2 ∈ (Ker (G))⊥ and GE (x1 ) = GE (x2 ). Then G (x1 − x2 ) = 0 and x1 − x2 ∈ Ker (G), which implies x1 = x2 . Thus GE is injective. The discussion above indicates that GE is an isomorphism. It 1 is not hard to prove the continuity of G− E , thus GE is also a homeomorphism. (ii) Since GE is an isomorphism, it is trivial that dim (Ker (G))⊥ = 3.

(

1 0

0 eJ2 t 0

⎞

0 0 ⎠, eJ3 t

t eiωt , eJ2 t = 1 0

)

(18)

(

0

e

iωt

)

, eJ3 t =

e−iωt 0

(

0

e

)

−iωt

.

(19)

We can calculate e−Aτ similarly. Substituting the expressions of eAt and e−Aτ into Eq. (6) (P0 = 0), we have r (T ) = G (f ) = ⟨f , φ1 ⟩L2 , ⟨f , φ2 ⟩L2 , ⟨f , φ3 ⟩L2

(

)⊤

/ω,

(20)

where φ1 = (sin (ω (T −τ )) , 2 (1 − cos (ω (T −τ ))) , 0)⊤

φ2 = (2 (cos (ω (T −τ ))− 1) , 4 sin (ω (T −τ ))− 3ω (T −τ ) , 0)⊤ .

(21)

φ3 = (0, 0, sin (ω (T −τ )))⊤ The following theorem indicates a basis of (Ker (G))⊥ . Theorem 3. The set {φ1 , φ2 , φ3 } defined by Eq. (21) is a basis of (Ker (G))⊥ . Proof. The condition G (f ) = 0 is equivalent to ⟨f , φi ⟩L2 = 0(i = 1, 2, 3), which implies φ1,2,3 ∈ (Ker (G))⊥ . Suppose there exist λ1 , λ2 , and λ3 ∈ R so that λ1 φ1 + λ2 φ2 + λ3 φ3 = 0. Then ( we have ) λ1 ⟨f , φ1 ⟩L2 + λ2 ⟨f , φ2 ⟩L2 + λ3 ⟨f , φ3 ⟩L2 = 0, ∀f ∈ L2 [0, T ] , R3 . Since G is surjective, we have λ1 x1 + λ2 x2 + λ3 x3 = 0, ∀ (x1 , x2 , x3 )⊤ ∈ R3 , which indicates λ1 = λ2 = λ3 = 0 and {φ1 , φ2 , φ3 } is linearly independent. Combining with dim (Ker (G))⊥ = 3, {φ1 , φ2 , φ3 } is a basis of (Ker (G))⊥ . To construct an orthonormal basis of (Ker (G))⊥ , we can utilize the Gram–Schmidt orthogonalization. Let

β1 = φ 1 ⟨φ2 , β1 ⟩L2 β1 , ⟨β1 , β1 ⟩L2 ⟨φ3 , β1 ⟩L2 ⟨φ3 , β2 ⟩L2 β3 = φ3 − β1 − β2 ⟨β1 , β1 ⟩L2 ⟨β2 , β2 ⟩L2 β2 = φ 2 −

(22)

and ei = βi /∥βi ∥L2 , (i = 1, 2, 3),

(23)

then {e1 , e2 , e3 } is an orthonormal basis of (Ker (G)) . Since Ker (G) is a Hilbert space, it also has} an orthonormal basis, { ( ) + 2 3 [0 ] which is denoted as en |n ∈ N , n ≥ 4 . As L , T , R = { } ⊥ + Ker (G ⊕ Ker G , the set e | n ∈ N is an orthonormal basis ) ( ( )) n ( ) of L2 [0, T ] , R3 . Since GE is an isomorphism between (Ker (G))⊥ and R3 , {GE (e1 ) , GE (e2 ) , GE (e3 )} is a basis of R3 . The explicit expressions of GE (ei ) , ei (i = 1, 2, 3) are given in the Appendix. ⊥

4

H. Sun and J. Li / Automatica 115 (2020) 108864

6. Structure of the reachable set

By Theorem 4, we obtain the following explicit expressions which describe the structure of the reachable set. Denote

The constructed orthonormal basis enables us to describe the structure of the reachable set. In this section, we demonstrate the theorem and explicit expressions which describe the reachable ( ) set. For f ∈ L2 [0, T ] , R3 , the Parseval equality (Yosida, 1978) indicates that

H11 H = (GE (e1 ) , GE (e2 ) , GE (e3 )) = H12 0

∥f ∥2L2 =

∞ ∑

⟨f , ei ⟩2L2 ,

(24)

(

√ Theorem 4. The reachable √ √ set under the constraint (i) (ii) 2E = R, or (iii) r ≤ 2E ≤ R is

0 0 H33

) ,

(32)

where H11 , H22 , H33 > 0 and H12 < 0 (see the Appendix). Then the reachable set can be expressed as

{

i=1

which associates the energy of f with its components along ei (i ∈ N+ ). We have the following theorem.

0 H22 0

} )⊤ ( (x, y, z )⊤ ∈ R3 |(x, y, z ) H −1 H −1 (x, y, z )⊤ ≤ R2 .

(33)

The boundary of the reachable set is

{

} )⊤ ( (x, y, z )⊤ ∈ R3 |(x, y, z ) H −1 H −1 (x, y, z )⊤ = R2 .

(34)

2E ≤ R, 7. Properties of the reachable set

{

} √ G (Ui ) = xGE (e1 )+ yGE (e2 )+ zGE (e3 ) | x2 + y2 + z 2 ≤ R , i = 1, 2, 3,

(25)

where U1,2,3 are defined by Eq. (10) and GE , e1,2,3 are defined in Sections 4 and 5. Proof. (i) i = 1 Since xGE (e1 ) + yG √ E (e2 ) + zGE (e3 ) = GE (xe1 + ye2 + ze3 ) and ∥xe1 + ye2 + ze3 ∥L2 = x2 + y2 + z 2 ≤ R, it is trivial that

{

√

xGE (e1 ) + yGE (e2 ) + zGE (e3 ) | x2 + y2 + z 2 ≤ R

} (26)

⊆ G (U1 ) .

Suppose f = f// + f⊥ ∈ U1 , where f// ∈ Ker (G) , f⊥ ∈ (Ker (G))⊥ . We have G (f ) = G (f⊥ ) = ⟨f , e1 ⟩L2 GE (e1 ) + ⟨f , e2 ⟩L2 GE (e2 ) + ⟨f , e3 ⟩L2 GE (e3 ) and

√ ⟨f , e1 ⟩2L2 + ⟨f , e2 ⟩2L2 + ⟨f , e3 ⟩2L2 ≤ ∥f ∥L2 ≤ R.

(27)

Thus we have

{

√

xGE (e1 ) + yGE (e2 ) + zGE (e3 ) | x2 + y2 + z 2 ≤ R

} (28)

⊇ G (U1 ) . According to the discussion above, we have

{

√

xGE (e1 ) + yGE (e2 ) + zGE (e3 ) | x2 + y2 + z 2 ≤ R

} (29)

= G (U1 ) . (ii) i = 2 Similar to (i), we have

{

√

xGE (e1 ) + yGE (e2 ) + zGE (e3 ) | x2 + y2 + z 2 ≤ R

} (30)

⊇ G (U2 ) .

√

Let f = xe1 + ye2 + ze3 + λe4 , where x2 + y2 + z 2 ≤ R and λ ∈ R. Then √ we have G (f ) = xGE (e1 ) + yGE (e2 ) + zGE (e3 ) and ∥f ∥L2 = x2 + y2 + z 2 + λ2 . We can select an appropriate λ ∈ R to have f ∈ U2 , thus

{

√

xGE (e1 ) + yGE (e2 ) + zGE (e3 ) | x2 + y2 + z 2 ≤ R

⊆ G (U2 ) .

} (31)

Combining the expressions (30) and (31), the proof is completed. (iii) i = 3 The proof is similar to that in (ii). Remark 1. By definition, it is trivial that G (U2 ) ⊆ G (U3 ) ⊆ G (U1 ). Theorem 4 indicates that G (U1 ) = G (U2 ) = G (U3 ). The reachable sets under the three different constraints are the same.

In this section, we use the explicit expressions obtained in Section 6 to study the shape, volume, surface area and other properties of the reachable set. The constraints on the control force are the same as those in Section 3. According to Remark 1, we can denote the reachable set as Ω = G (U1 ) = G (U2 ) = G (U3 ). (i) Denote B (0, R) as the ball in R3 , where the radius R is given by the energy constraint in Section 3. We define the map H : B (0, R) → Ω , (x, y, z )⊤ ↦ → H (x, y, z )⊤ . The map H is a diffeomorphism, which enables us to calculate the volume and surface area of the reachable set. By the substitution of variables, the volume of the reachable set is V (Ω ) =

2π

∫ 4

R

∫

|det (H )| r 2 sin θ drdθ dφ 0

0

0

=

π

∫

.

(35)

π R3 H11 H22 H33

3

⊤ Let A ) = H (R sin θ cos φ, R sin θ sin φ, R cos θ) and ∇ A = (∂θ A, ∂φ A , then the surface area of the reachable set is

S (∂ Ω ) =

2π

∫ 0

∫

π

√

det (∇ A)⊤ ∇ A dθ dφ.

(

)

(36)

0

(ii) By the expressions of GE (e1 ) , GE (e2 ) , GE (e3 ) in the Appendix, we can study the shape of the reachable set. The directions of GE (e2 ) , GE (e3 ) are the same as those of y-axis, z-axis respectively. By Eqs. (A.2) and (A.4), the components of GE (e1 ) along x-axis, y-axis, and z-axis are H11 > 0, H12 < 0, H13 = 0 respectively. Suppose xGE (e1 )+ yGE (e2 )+ zGE (e3 ) is a point in the reachable set. By Theorem 4, xGE (e1 ) + yGE (e2 ) − zGE (e3 ) is also a point in the reachable set. Thus the reachable set is symmetric with respect to the xy-plane. The direction of GE (e1 ) deviates from x-axis by an acute angle. By Theorem 4 again, the shape of the reachable set is an ellipsoid with one deviated axis. Define ϕ = arctan (−H12 /H11 ), which is the angle between GE (e1 ) and x-axis. By direct calculation, we have

ϕ = arctan

(

12 (sin (T ω) − T ω)2 3 sin (2T ω) − 32 sin (T ω) + 26T ω

)

.

(37)

The following theorem describes the asymptotic behavior of the angle ϕ when the control time T changes. The proof of Theorem 5 only utilizes elementary methods in calculus, and thus it is omitted. Theorem 5. Suppose the angular velocity ω is fixed, then (i) ϕ is a monotonically increasing function of T , (ii) limT →0+ ϕ (T ) = 0, (iii) limT →+∞ ϕ (T ) = π/2.

H. Sun and J. Li / Automatica 115 (2020) 108864

5

Remark 2. Note that Theorem 5(iii) is only a theoretical result. Since Hill’s equations describe the linearized relative motion, the reachable set does not coincide with the practical one when T is too large. 8. Numerical verification In this section, we use numerical methods to verify the theoretical conclusions obtained. √ (i) We calculate the reachable set under the constraint 2E ≤ R by the numerical method and Theorem 4 respectively. We compare the results } to do the verification. The set {1, sin (nx) , cos (nx) |n ∈ N+ is an orthogonal basis of L2 ([−π, π ] , R) (Yosida, 1978). By translation and stretching transformation, the set

{√

2 T

⋃

(

sin nπ

(

2t T

} )) √ ( ( )) 2 2t + −1 , cos nπ − 1 |n ∈ N T

Fig. 2. Let R = 0.0361 m s−3/2 , ω = 9.2002 × 10−4 rad/s, N = 6 and T = 1300 s. The 104 red points are the final positions under randomly selected control forces, and the blue surface is the boundary of reachable set calculated by Theorem 4. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

T

{√ } 1

T (38) is an orthonormal basis of L ([0, T ] , R), and the control force can be approximated by a combination of finite elements in the basis. Suppose the control force f (t ) = (f1 (t ) , f2 (t ) , f3 (t ))⊤ is 2

√ fi ( t ) =

1 T

√

a(i) +

N 2∑

T

√

N 2∑

T (i)

(

bn(i) sin nπ

(

n=1

(

cn cos nπ

n=1

(

2t T

2t T

)) −1

+

))

,

(39)

−1

where i = 1, 2, 3 and N is a given positive integer. The control force is required to satisfy ) ( N ( 3 ∑ ( (i) )2 ∑ ( (i) )2 ( (i) )2 ) 2 ∥f ∥L2 = a + bn + cn ≤ R2 . (40) i=1

n=1

We randomly select 104 groups of coefficients a(i) , bn(i) , cn(i) |i = 1, 2, 3, 1 ≤ n ≤ N

{

}

(41)

which satisfy the inequality (40), and calculate the corresponding final positions by Eq. (20). As shown in Fig. 2, all the final positions under randomly selected control forces are in the reachable set calculated by Theorem 4. The numerical results verify the correctness of theoretical conclusions. (ii) The reachable sets with different control time are calculated by Theorem 4, and the results are shown in Figs. 3, 4, and 5. The values of R, ω are the same as Fig. 2. The numerical results indicate that, as T increases, the angle between GE (e1 ) and x-axis increases and the reachable set deviates more from x-axis, which is consistent with Theorem 5. 9. Conclusion Using the methods in Hilbert space theory, we give Theorem 4 and the explicit expressions (33), (34) for computing the reachable set for spacecraft relative motion under different constraints √ of energy. The reachable sets √ √ under the constraint (i) 2E ≤ R, (ii) 2E = R, or (iii) r ≤ 2E ≤ R are the same. The shape of the reachable set is an ellipsoid with one deviated axis. We also calculate the volume and the surface area of the reachable set. The numerical methods are used to verify the conclusions above. This article thoroughly solves the computation of the reachable set for spacecraft relative motion under the energy-limited

Fig. 3. Projection of the reachable set to xy-plane when T = 500 s.

constraint. The explicit expressions obtained in this article can be used to calculate the accurate reachable set without any computational loads. This enables the spacecraft to maintain high-level and immediate awareness about its operational capabilities and limitation. Acknowledgments We gratefully acknowledge the valuable cooperation of Yi Zhu and Meirong Zhang (Zhou Pei-Yuan Center for Applied Mathematics, Tsinghua University). Appendix. Calculation and explicit expressions We calculate the explicit expressions of GE (ei ) , ei (i=1, 2, 3). Note that ei , φi , βi (i=1, 2, 3) have been given by Eqs. (21), (22), and (23). Denote η = ∥β2 ∥L2 and

ξ=

⟨φ2 , β1 ⟩L2 12 (sin (T ω)− T ω)2 =− , ⟨β1 , β1 ⟩L2 3 sin (2T ω)− 32 sin (T ω)+ 26T ω

(A.1)

6

H. Sun and J. Li / Automatica 115 (2020) 108864

Fig. 4. Projection of the reachable set to xy-plane when T = 1500 s.

Fig. 5. Projection of the reachable set to xy-plane when T = 2500 s.

ei = (ei1 , ei2 , ei3 )⊤ , GE (ei ) = (Hi1 , Hi2 , Hi3 )⊤ , i = 1, 2, 3.

(A.2)

H23 = H31 = H32 = 0

By direct calculation, we have

) ( ) 1 ( η = 8ξ (3 − 4ξ ) sin2 (T ω) + 3 ξ 2 − 4 sin (2T ω) − 4ω ( ) 32 T 38 + 3T 2 ω2 − 12ξ sin (T ω) + sin (T ω) + ( T

6ξ T ω +

13 2

ξ − 48 sin 2

2

(

Tω

ω ))

3 sin (2T ω) − 32 sin (T ω) + 26T ω

ω

e13 = 0

(2 (cos (ω (T − τ )) − 1) − ξ sin (ω (T − τ )))

(4 sin (ω (T − τ )) − 3ω (T − τ )) − η 2ξ (1 − cos (ω (T − τ ))) η

1

e23 = e31 = e32 = 0 e33 = 2 sin (ω (T − τ ))

ω

√

2T ω − sin (2T ω)

1 √

ω (3 sin (2T ω) − 32 sin (T ω) + 26T ω) )2 sin (T ω) H12 = − 6 T − · ω √ ω 3 sin (2T ω) − 32 sin (T ω) + 26T ω H13 = H21 = 0 H11 =

2ω2

(

H22 =

1

ω2 η (3 sin (2T ω) − 32 sin (T ω) + 26T ω)

)

)

2

2

)

4

(A.5)

4

References

3 sin (2T ω) − 32 sin (T ω) + 26T ω

e22 =

(

16T ω sin (T ω) 3T ω − 71 + 42T ω

e12 =4 (1 − cos (ω (T − τ ))) ·

1

ˆ 22 =2 cos (2T ω) 54T 2 ω2 − 247 + 624T 2 ω2 cos (T ω)+ H

It is clear that H11 , H33 > 0. Now we indicate H22 > 0 ˆ 22 > 0. Let x = T ω and the and it is sufficient to prove H ˆ iterm in H22 (x) with the highest order is 42x4 . Thus we have ˆ 22 (x) → +∞ (x → +∞). We also have Hˆ 22 (x) ̸= 0, ∀x > 0, H otherwise GE (e2 ) = 0. By the intermediate value theorem, we ˆ 22 (x) > 0, ∀x > 0. have H

ω

η

(A.4)

where (A.3)

(

2

e21 =

ω (2T ω − sin (2T ω)),

(

e11 =2 sin (ω (T − τ )) ·

√

ω2

3T ω sin (2T ω) 3T 2 ω2 − 140 + 494 + 256T 2 ω2 .

and

√

1 √

H33 =

2

ˆ 22 H

Blanchini, F. (1999). Survey paper: Set invariance in control. Pergamon Press. Chen, Xin, Erika, & Sriram (2013). Flow*: An analyzer for non-linear hybrid systems. In International conference on computer aided verification (pp. 258–263). Fehse, W. (2003). Automated rendezvous and docking of spacecraft: Space and ground system setup. Cambridge University Press. Frehse, G., Guernic, C. L., et al. (2011). SpaceEx: Scalable verification of hybrid systems. In International conference on computer aided verification (pp. 379–395). Hill, G. W. (1878). Researches in the Lunar theory. American Journal of Mathematics, 1, 5–26. Kurzhanski, A. B, & Varaiya, P. (2000). Ellipsoidal techniques for reachability analysis. In International workshop on hybrid systems: Computation and control (pp. 202–214). Lee, S., & Hwang, I. (2018). Reachable set computation for spacecraft relative motion with energy-limited low-thrust. Aerospace Science and Technology, 77, 180–188. Schlanbusch, R., Kristiansen, R., & Nicklasson, P. J. (2011). Spacecraft formation reconfiguration with collision avoidance. Automatica, 47(7), 1443–1449. Wen, C., & Gurfil, P. (2016). Relative reachable domain for spacecraft with initial state uncertainties. Journal of Guidance, Control and Dynamics, 39(3), 1–12. Xue, Dan, Li, Junfeng, & Baoyin, Hexi (2015). Reachable domain for spacecraft with a single impulse. Journal of Guidance, Control and Dynamics, 33(3), 934–942. Yosida, K. (1978). Functional analysis. Springer-Verlag.

H. Sun and J. Li / Automatica 115 (2020) 108864 Hao Sun received the bachelor’s degree in mathematics and physics from Tsinghua University, China, in 2018. After graduation, he is a Ph.D. candidate in Yau Mathematical Sciences Center, Tsinghua University. He majors in fundamental mathematics and his research interests include dynamical systems, celestial mechanics, and geometric analysis.

7 Junfeng Li was born in Heilongjiang, China. He received the bachelor’s degree in mechanics from Peking University, China, in 1987, and the Ph.D. in mechanics from Moscow State University, Russia, in 1993. He worked at the Surrey space center, UK, from 1998 to 1999. Since 1999, he works as a professor at Tsinghua University in Beijing, China. His research interests include stability of motion, astrodynamics, multi-body dynamics, and liquid sloshing.