- Email: [email protected]
Modeling and analysis of nonlinear spacecraft relative motion via harmonic balance and Lyapunov function
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.1 (1-11)
Aerospace Science and Technology ••• (••••) ••••••
1
Contents lists available at ScienceDirect
67 68
2 3
Aerospace Science and Technology
4
69 70 71
5
72
6
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13
Modeling and analysis of nonlinear spacecraft relative motion via harmonic balance and Lyapunov function ✩
16 17
81
Mechanical and Aerospace Engineering Department, Naval Postgraduate School, Monterey, CA, United States of America
83
82 84 85
19
a r t i c l e
i n f o
a b s t r a c t
86 87
21 22 23 24 25 26 27 28 29 30 31 32
79
Ayansola D. Ogundele
18 20
78 80
14 15
77
Article history: Received 14 October 2019 Received in revised form 27 October 2019 Accepted 1 February 2020 Available online xxxx Communicated by Efstathios Bakolas Keywords: Relative motion Cubic approximation Harmonic balance Lyapunov function Total energy
The harmonic balance method, an approximate method for the study of nonlinear oscillating systems described by ordinary nonlinear differential equations, is applied to the cubic polynomial approximation of the relative motion to develop approximate periodic solutions of the motion. Firstly, the cubic nonlinear terms neglected in the formulation of linear time invariant Clohessy-Wilshire (CW) equations are obtained from the original nonlinear equation of motion. Secondly, using the CW solutions as the generating solutions, harmonic balance technique is applied to the cubic nonlinear terms and correction terms are obtained. The new and improved harmonically linearized model contains CW equations and the correction terms. Using total energy as a Lyapunov function, the stability of the new model is determined. The additional terms in the new harmonically linearized spacecraft relative motion enabled us to obtain excellent approximate periodic solutions of the nonlinear equations as evidenced from the numerical simulation results. © 2020 Elsevier Masson SAS. All rights reserved.
88 89 90 91 92 93 94 95 96 97 98
33
99
34
100
35
101
36
102
37 38
1. Introduction
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
George Hill (1878) in the nineteenth century discovered that a great insight into the dynamics of nearby objects can be gained by linearizing the relative equations of motion and ignoring the longterm effects of the natural perturbations in the context of celestial mechanics [1]. In 1960, W. R. Clohessy and R. S. Wiltshire used this approach when researching on terminal guidance for satellite rendezvous in the context of astrodynamics [2]. The relative motion dynamics between two spacecraft, for example Chief and Deputy spacecraft, are generally analyzed using Clohessy-Wilshire (HCW) equations. In this model, a reference frame centered on a reference spacecraft, called Chief spacecraft, that orbits the Earth in a circular trajectory is defined. The assessment of the operational capability of the spacecraft in formation flying can be done using the knowledge on a reachable set for their relative motion [3]. Proper modeling of relative motion dynamics is desirable as this will enhance spacecraft autonomy which is required by the next generation of spacecraft [4,5]. It was established that [6] the relative position dynamics is more suitable for actual measurement circumstances and applicable for arbitrary orbital forms.
59 60 61 62 63 64 65 66
✩
This document is a revised version of the paper presented at 2016 AAS/AIAA Spaceflight Mechanics Meeting. E-mail address: [email protected]. https://doi.org/10.1016/j.ast.2020.105761 1270-9638/© 2020 Elsevier Masson SAS. All rights reserved.
The motion of Deputy spacecraft moving relative to the Chief is studied from a reference frame fixed at the center of the reference spacecraft called Local Vertical Local Horizontal frame [2,7–12]. The orbital radius of the reference spacecraft is assumed to be much larger than the relative separation between the spacecraft. Over a long time this approach has significant errors. Because of this limitation there have being some modifications in one way or the other over a period of time. Battin uses universal variables which allows it to be used for elliptical and hyperbolic orbits [13,14]. Also, Carter and Humi modified the dynamics for eccentric reference orbits [15]. London [16] developed and solved second order polynomial approximation of satellite equations of motion but a better accuracy can be obtained using higher order polynomial. In [17] detailed modeling and analysis of spacecraft relative motion dynamics actuated by inter-satellite non-contacting force is given. Survey of the technological difficulties and developmental perspectives in spacecraft dynamic analysis are highlighted in [18]. A comprehensive survey and assessment of spacecraft relative motion can be found in [19]. The harmonic balance (HB) method, or harmonic linearization, is a powerful technique for the analysis of high-frequency nonlinear oscillators and for obtaining analytic solutions to nonlinear ordinary differential equations describing oscillatory systems [20,21]. The method became prominent in 1990s because of its ability to calculate and capture the steady-state behavior of nonlinear systems directly and high applicability. Harmonic balance
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.2 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
method is based on trigonometric functions and it is possible to formulate this procedure using any complete set of periodic functions such as Jacobi elliptic functions. It can be adapted to construct a linear substitute for the original equation. The difficulties in using this method are that calculating the amplitudes and the angular frequency may become algebraically intensive and existing formulations of the harmonic balance procedure do not allow it to be applied to non-conservative oscillators because these types of solutions have solutions involving transient behaviors. In order to obtain a consistent solution using the method one needs to either know a great deal about the solution a priori or to carry enough terms in the solution and check the order of the coefficients of all the neglected harmonics, otherwise one might obtain inaccurate approximation [22,23]. Before applying harmonic balance the governing equations are firstly linearized to a lower order polynomial so as to prevent solving complicated nonlinear algebraic equations. This leads to linear algebraic equations instead of nonlinear one. The harmonic linearization has been applied successfully to several nonlinear systems [24–26]. The main contributions of this paper are itemized as follows:
22 23
• A new harmonically linearized spacecraft relative motion is
24 25
•
26 27 28
•
29 30 31
•
32 33 34
•
developed making use of assumed generating solution. The method is amenable based on the fact that it can be applied to solve differential equations whose nonlinear terms are not small. The resulted harmonic balance linearized form gave a better approximation of periodic solutions of the original nonlinear term. The new model is useful for analyzing spacecraft formation flying, rendezvous and proximity operations. Using total energy as a Lyapunov function the developed harmonic balance model stability was investigated.
35 36 37 38 39 40 41 42 43 44 45 46 47 48
In this paper, we developed harmonically linearized spacecraft relative motion using harmonic balance technique. To the best of our knowledge, this method is rarely used in the papers dealing with the spacecraft relative motion dynamics. The approach is a two stage procedure. The first approach involves linearization of the spacecraft nonlinear equation of motion into a cubic polynomial model using Binomial theorem and the second approach involves the application of harmonic linearization technique into the cubic model. This leads to the formulation of new harmonic balance relative equation of motion consisting of improved, three linear, second order ODEs. We investigated the stability of the model using total energy as Lyapunov function and it was found to be stable.
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
2. Review of harmonic balance method The review of harmonic balance method is adapted from the References [27,28,22,20,21,29]. This method is most commonly used to study the nonlinear vibration problems. It may be used to determine the approximate periodic and quasi-periodic oscillations, periodic and quasi-periodic conditions in automatic control theory, as well as stationary conditions, and in the studies of their stability. If a periodic solution does exist, it may be sought in the form of a Fourier series, whose coefficients are determined by requiring the series to satisfy the equation of motion. Consider a nonlinear perturbed oscillator differential equation
x¨ + ω x + ε f (x, x˙ ) = 0, 2
(0 < ε 1)
(1)
where ε is a small parameter. Harmonic linearization technique involves the replacement of the nonlinear forcing function F (x, x˙ ) = ε f (x, x˙ ) by the equivalent linear force [27]
67 68 69 70 71 72 73 74 75 76 77 78 79
Fig. 1. Steps followed in finding an approximate model for the nonlinear equation using Harmonic Balance method.
80 81 82
F l (x, x˙ ) = kx + λ˙x
(2)
83 84
The parameters k and λ can be computed from
ε ∫ f (a cos ψ, −aω sin ψ) cos ψ dψ πa 0 ε 2π ∫ f (a cos sψ, −aω sin ψ) sin ψ dψ λ (a) = − π aω 0
k (a) =
85 86
2π
87 88
(3)
where ψ = ωt + θ . If x = acos (ωt + θ), a= constant, ω = constant, the nonlinear force F (x, x˙ ) is a periodic function of time, and its Fourier series expansion contains an infinite number of harmonics, having the frequencies nω , n = 1, 2, . . .. The forcing function F (x, x˙ ) has the form
F (x, x˙ ) =
∞
89 90 91 92 93 94 95 96 97
F n cos (nωt + θn )
(4)
n =0
98 99 100
The term F 1 cos(ωt + θ1 ) is called the fundamental harmonic of the expansion. The amplitude and the phase of the linear function F l coincide with the respective characteristics of the fundamental harmonic of the nonlinear force. For the differential equation
101 102 103 104 105
x¨ + ω2 x + F (x, x˙ ) = 0
(5)
106 107
which is typical in the theory of quasi-linear oscillations, the harmonic balance method consists in replacing F (x, x˙ ) by the equivalent linear force to give
x¨ + λ˙x + k1 x = 0
108 109 110 111
(6)
112
where k1 = k + ω2 , λ is the equivalent damping coefficient and k1 is the equivalent elasticity coefficient. In the harmonic balance method the frequency of the oscillation depends on the amplitude (through the quantities k and λ). By traditional linearization about the origin we have
113
x¨ + ω x = 0
119
2
(7)
Comparing the harmonic linearization result in Eq. (6) with Eq. (7) it can be easily seen that harmonic linearized equation has correction terms. In this work, harmonic balance method is adapted to construct a new and improved linearized model of the spacecraft relative motion. Fig. 1 illustrates the steps followed in finding an approximate model for a nonlinear system. In developing the harmonicbalance model for spacecraft relative motion, polynomial approximation is applied to the nonlinear equations in LVLH components and it produced cubic equations of motion containing linear and nonlinear terms. The linear equations correspond to the ClohessyWiltshire equations. For the harmonic linearization, the nonlinear
114 115 116 117 118 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.3 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
3
The LHS of Eq. (10) shows the differentiation in inertial frame while that in the RHS shows the differentiation in Hill frame. Therefore, differentiating Eq. (8) twice and using Eq. (10) gives the kinematics equation of motion as
1 2 3 4
¨ + ρ¨ + 2ω × ρ˙ + ω × (ω × ρ ) r¨ = R
6 7 8
(11)
For the case with the chief in circular orbit the mean motion is ˙f = n = μ/ R 3 . This reduces Eq. (11) to
9 10
11
r¨ = x¨ − 2ny − n2 x −
12
μ ˆ
i + y¨ + 2nx˙ − n2 y ˆj + z¨ kˆ
R2
(12)
Using Eq. (12) in the second part of Eq. (9) gives the second order nonlinear differential equations of relative motion as
14 15 16
x¨ − 2ny − n2 x −
17 18
μ R2
y¨ + 2nx˙ − n y = −
20
Fig. 2. Relative Motion of Deputy spacecraft with respect to Chief Spacecraft.
z¨ = −
22
25 26
terms are evaluated along the generating solution to produce linear correction terms. The correction terms are combined with the original linear equations (CW equations) to produce the harmonic balance model of the spacecraft relative motion.
27 28
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
˜
58 59 60 61 62 63 64 65 66
72 73 74 75 76 77 78
horizontal frame (xyz) with unit vectors
ˆi, ˆj, kˆ , referred to as
Hill frame, is attached to the chief and it is rotating with it. The x axis is aligned along the outward radial direction, the z axis is aligned with the direction perpendicular to the orbit plane and it is parallel to the orbit momentum vector, and the y axis is normal to the other two axes and completes the right-hand coordinate system. The in-plane motion is defined by x and y while the out-of-plane motion is defined by the z axis. Expressed in the Hill frame, the deputy’s position relative to the chief is denoted as ρ = xˆi + yˆj + zkˆ = (x, y , z)T . The unperturbed rotating Hill (LVLH) frame relative to the inertial frame has constant angular velocity ˙ = ¨f kˆ = 0. From of ω = ˙f kˆ = Rh2 kˆ and angular acceleration of ω Fig. 1, R is the chief spacecraft position vector relative to the Earth and can be expressed as R = R ˆi in the Hill frame. The angular mo˙ = R 2 ˙f kˆ and its mentum expressed in the Hill frame is h = R × R 2˙ magnitude is h = R f . The deputy spacecraft position vector r relative to the center of the gravitational field is
r=R+ρ
r3
˜
(8)
˜
z + z˜ = V x˜ , y˜ , z˜
(13)
88
¨ =− R r¨ = −
R3
(14)
μ r3
r
(9)
2 Nd ρ
dt 2
= ρ¨ + 2ω × ρ˙ + ω × (ω × ρ )
(10)
94 95 97
99
2
3 2
y˜ + 2
3
100
z˜ − 6x˜ y˜ − 6x˜ z˜ − 4x˜ 2
2
2 3 2 3 3 2 U x˜ , y˜ , z˜ = 3 y˜ x˜ − 6 y˜ x˜ + y˜ z˜ − y˜ 2 2 3 3 V x˜ , y˜ , z˜ = 3 z˜ x˜ − 6 z˜ x˜ 2 + z˜ y˜ 2 − z˜ 3 2 2
2
3
101 102
(15)
105 106 107 108 109 110 111 112
x˜ − 2 y˜ − 3x˜ = 0 y˜ + 2x˜ = 0
103 104
Eq. (14) is the nonlinear cubic approximation model of spacecraft relative motion with chief in circular orbit. For ε = 0, that is neglecting the nonlinear terms, the system in Eqs. (14) reduced to the linear equations
113
(16)
114 115
˜
z + z˜ = 0
116
These equations are the well-known normalized Clohessy-Wilshire (CW) Equations. They are coupled equations of motion in x and y. Using linear theory, the homogeneous equations have the solutions [9,29]
118
117 119 120 121 122 123 124
(17)
127 128
τ = nt and
c 1 = 2x˜ 0 + y˜ 0 , c 2 = x˜ 0 , c 3 = −3x˜ 0 − 2 y˜ c 4 = −2x˜ 0 + y˜ 0 , c 5 = z˜ 0 , c 6 = z˜ 0
125 126
z˜ = c 5 sin τ + c 6 cos τ where,
Using transport theorem we obtain
91
96
y˜ = −3c 1 τ + 2c 2 cos τ − 2c 3 sin τ + c 4
R
90
93
x˜ = 2c 1 + c 2 sin τ + c 3 cos τ
μ
89
92
The Keplerian acceleration of the chief and deputy spacecraft are
86
98
T x˜ , y˜ , z˜ = −3x˜ +
85 87
z
where,
81 83
y
r3
μ
80
84
μ
y + 2x = U x˜ , y˜ , z˜
In this section, the cubic polynomial approximation model of deputy spacecraft motion relative to the chief spacecraft in circular orbit is derived. The derivation follows the approach in [10,29,28]. Fig. 2 shows relative motion of deputy spacecraft with respect to the chief spacecraft in circular orbit. A local-vertical/local
( R + x)
r3
x˜ − 2 y˜ − 3x˜ = T x˜ , y˜ , z˜
3. Nonlinear cubic approximation model of relative motion
56 57
70
82
μ
Applying Binomial theorem expansion to Eq. (13) and normalizing the resulting equations we have
29 30
=−
2
19
24
69
79
13
23
68
71
5
21
67
129 130
(18)
131 132
JID:AESCTE
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.4 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
4
1 2
4. Harmonic linearization of nonlinear cubic approximation model
Substituting the linear approximations in Eqs. (23) and (24) into Eq. (23) and use the results in the radial part of Eq. (14) gives
3 4 5 6 7 8 9
In this section, harmonic linearization is applied to the cubic approximation model of spacecraft relative motion using the bounded, steady-state solutions of CW equations as the generating solutions. In order to obtain bounded solutions, the secular term in the along-track is eliminated by setting c 1 = 0 and this leads to
10
x˜ = c 2 sin τ + c 3 cos τ
11
y˜ = 2c 2 cos τ − 2c 3 sin τ + c 4
12 13
z˜ = c 5 sin τ + c 6 cos τ
14 15 16 18 19 20 21 22 23
The radial nonlinear part, from Eq. (14), is given by
T x˜ , y˜ , z˜ = −3x˜ +
cos2 τ =
26 27 28
2
y˜ + 2
3
z˜ − 6x˜ y˜ − 6x˜ z˜ − 4x˜ 2
2
2
2
3
cos
3
τ=
1 2 1 4
1
sin2 τ =
(1 + cos 2τ ),
2 3
(3 cos τ + cos 3τ ),
sin
τ=
1
(3 sin τ − sin 3τ )
4
31 32
(21)
33 34 35
The nonlinear terms are linearized as follows
x˜ = 2c 2 c 3 sin τ cos τ 2
+ c 22 sin2
τ
+ c 32 cos2
τ
Substituting Eq. (21) into Eq. (22) and eliminating the higher order harmonic of time yields
41
Similar simplification is performed for the remaining terms as
38 39
42 43 44 45 46 47 48
y˜ = 2
τ
51 52 53 54 55
≈
56
59 60 61 62 63 64
x˜ z˜
2
c 22
τ
τ
τ
τ
τ
τ
τ
τ
+ c 32
+ c 42
= c 3 c 62 cos3 +c 3 c 52 cos
τ
τ
U x˜ , y˜ , z˜ = 3 y˜ x˜ − 6 y˜ x˜ 2 +
3 2
y˜ z˜ 2 −
3
τ sin τ
+2c 3 c 5 c 6 cos τ sin τ 1 ≈ (2c 2 c 5 + 2c 3 c 6 ) z˜ + c 52 + c 62 x˜
4 x˜ 3 = c 33 cos3 τ + c 23 sin3 τ + 3c 2 c 32 cos2 τ sin τ
3 2 ≈ c 2 + c 32 x˜ 4
80 81 82 83 84 86
2
88 89
(27)
y˜ x˜ = 4c 2 c 3 sin τ cos
2
2
3
2
95 97
100 101
3
τ + 2c2 sin τ cos τ + 2c2 c3 cos τ
102 103 104
−2c 3 sin τ cos τ + 2c 2 c 3 c 4 sin τ cos τ
105
+c 2 2 c 4 sin2 τ + c 3 2 c 4 cos2 τ 1 2 = c 2 c 4 + c 3 2 c 4 + c 2 2 + c 3 2 y˜ 2
94
99
2
3
93
98
The remaining nonlinear terms are harmonically linearized as 2
92
96
(28)
−2c 3 sin τ cos τ + c 2 c 4 sin τ + c 4 c 3 cos τ ≈ c 4 x˜ 2
90 91
2
106 107 108 109 110
3
+2c 2 c 6 cos τ − 2c 3 c 5 sin τ
111
−4c 3 c 5 c 6 sin2 τ cos τ − 2c 3 c 6 2 sin τ cos2 τ
113
2
+c 4 c 6 2 cos2 τ 1 2 = c 5 + c 6 2 c 4 + c 5 2 + c 6 2 y˜ 4
+4 (c 2 c 6 − c 3 c 5 ) z˜ y˜ 3 = 2c 2 c 4 2 cos τ + 8c 2 2 c 4 cos2 τ + 2c 2 c 4 2 cos τ 2
+8c 2 c 4 cos τ − 8c 2 c 3 c 4 sin τ cos τ
2
79
y˜ 3
+c 4 c 5 sin τ + 2c 4 c 5 c 6 sin τ cos τ
(24)
77 78
Harmonic linearization of Eq. (27) is performed as follows. Substituting the trigonometric identities in Eq. (21) into Eq. (27), eliminating the higher order harmonics and rearrange the result in the form of assumed linear solution yields
2
τ
τ sin τ + 2c2 c5 c6 cos τ sin2 τ
+3c 22 c 3 cos τ sin2 τ
+ 3c 62
76
87
2
+ c 2 c 62 cos2
+ 3c 52
The along-track nonlinear part, from Eq. (14), is
2
τ
x˜
+ c 2 c 52 sin3 2
+ 18c 42
4 y˜ z˜ 2 = 2c 2 c 5 2 cos τ sin2 τ + 4c 2 c 5 c 6 sin τ cos2 τ
τ
τ
65 66
(23)
2 + c 3 c 42 cos τ
+c 2 c 42 sin + 4c 22 c 3 cos3 + 4c 2 c 32 sin3 +4c 22 c 4 cos sin − 4c 32 c 4 cos sin −8c 2 c 32 cos2 sin − 8c 22 c 3 cos sin2 +4c 2 c 3 c 4 cos2 − 4c 2 c 3 c 4 sin2
+ 6c 32
3
τ − 4c3 c4 sin τ + 4c2 c4 cos τ
x˜ y˜ 2 = 4c 23 cos2 τ sin τ + 4c 33 cos τ sin2 τ
50
58
τ
+ 4c 32 sin2
4
6c 22
(26)
−4c 2 c 3 2 sin2 τ cos τ − 2c 2 2 c 3 sin3 τ
−8c 2 c 3 cos τ sin τ + c 42 ≈ −2c 4 y˜ + 2c 22 + 2c 32 + 3c 42 1 2 z˜ 2 = c 62 cos2 τ + c 52 sin2 τ + 2c 5 c 6 cos τ sin τ ≈ c 5 + c 62
49
57
4c 22 cos2
75
85
(22)
40
37
1
74
y˜ x˜ = 2c 2 2 sin τ cos τ + 2c 2 c 3 cos2 τ − 2c 2 c 3 sin2 τ
1 2 1 2 x˜ 2 = c 2 c 3 sin 2τ + c 3 − c 22 cos 2τ + c 2 + c 32 2 2 1 2 2 ≈ c + c3 2 2
36
73
2 a42 = 3c 4 ,
b4 =
70 72
4.2. Harmonic linearization of along-track equation of motion
(1 − cos 2τ ),
29 30
−18c 22 − 18c 32 − 12c 42 − 3c 52 − 3c 62 ,
(20)
Harmonic linearization of each of the nonlinear terms is performed by substituting the generating solution in Eq. (19) into the nonlinear terms and using the following trigonometric identities:
24 25
3
2
1
Eq. (25) is the radial harmonically linearized relative motion. The additional terms in x, y , z and the nonhomogeneous term b4 are the radial corrections to the equation of motion. They are approximations of the cubic powers that were retained in the expansion of the nonlinear equations of motion.
4.1. Harmonic linearization of radial equation of motion
17
(25)
71
where,
a43 = −3 (c 2 c 5 + c 3 c 6 ), (19)
68 69
x˜ − 2 y˜ − (3 + a41 ) x˜ − a42 y˜ − a43 z˜ = b4
a41 =
67
−16c2 2 c 3 sin τ cos2 τ + 8c 2 3 cos3 τ 2
2
2
+8c 2 c 3 sin τ cos τ − 2c 3 c 4 sin τ −8c 2 c 3 c 4 sin τ cos τ + 8c 3 2 c 4 sin2 τ 2
2
+16c2 c 3 sin τ cos τ = 3c 2 2 c 4 + 3c 3 2 c 4 − 2c 4 3 + 3c 2 2 + 3c 3 2 + 3c 4 2 y˜
112
(29)
114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.5 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
1 2 3 4 5 6 7 8 9 10
Substituting these linear approximations into Eq. (27) and using the results in the along-track part of Eq. (14) leads to the harmonically linearized along-track equation of motion
The additional terms in x, y , z and the nonhomogeneous term b3 (c ) are the corrections to the cross-track direction.
y˜ + 2x˜ − a51 x˜ − a52 y˜ − a53 z˜ = b5
4.4. Model summary
13 14 15 16 17 18 19
(30)
where,
3 a52 = −16c2 2 − 16c3 2 − 12c4 2 + c 5 2 + c 6 2 , a53 = b5 =
8 3
(31)
8
−16c4 c 2 2 − 16c4 c 3 2 + c 4 c 5 2 + c 4 c 6 2 + 8c 4 3
The additional terms in x, y , z and the nonhomogeneous term b5 (c ) are the corrections to the along-track direction. 4.3. Harmonic linearization of cross-track equation of motion
20 21
The cross-track nonlinear part, from Eq. (14), is
22 23 24 25 26 27 28 29 30 31
V x˜ , y˜ , z˜ = 3 z˜ x˜ − 6 z˜ x˜ 2 +
2
z˜ y˜ 2 −
3 2
z˜ 3
(32)
2
2
35
3
4 z˜ y˜ 2 = 4c 5 c 2 2 cos2 τ sin τ + 4c 6 c 2 2 cos3 τ
38
(33) 2
39
+4c 2 c 4 c 5 cos τ sin τ + 4c 2 c 4 c 6 cos τ
40
−4c 3 c 4 c 5 sin2 τ − 4c 3 c 4 c 6 sin τ cos τ
41
2
44
2
47 48
z˜ = c 5 sin 3
3
52
55
2
2
z − a61 x˜ − a62 y˜ + (1 − a63 ) z˜ = b6
57
where,
2
(34)
63 64 65 66
82
⎡
⎤ ⎡ 0 0 0 1 0 0 0 0 0 0 0 1 0⎥ ⎢ 0 ⎢ ⎢ ⎢ 0 0 0 0 0 1⎥ 0 ⎢ ⎥ ⎢ x˜ = ⎢ + 3 + a41 (c) a42 (c) a43 (c) 0 2 0 ⎥ ⎢ b 4 ( c) ⎣ ⎦ ⎣ a51 (c) a52 (c) a53 (c) b 5 ( c) −2 0 0 a61 (c) a62 (c) a63 (c) − 1 0 0 0 b 6 ( c) x˜ = A (c) x˜ + B (c)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦
a62 =
3
(c 2 c 6 − c 3 c 5 ) 2 1 a63 = 12c 4 2 − 9c 5 2 − 9c 6 2 , 8 1 12c 4 2 − 9c 5 2 − 9c 6 2 b6 = 8
83 84 85 86 87 88 89 90 92 93 94
In this section, total energy is applied as a Lyapunov function to investigate the stability of newly developed harmonic balance model of spacecraft relative motion. Several authors have studied vibrations and stability of gyroscopic conservative systems [30–32]. General discussions and developments on gyroscopic systems, and a variety of investigations ranging from important practical problems to new general theories can be found in [33–36]. In 1989, Koditschek applied total energy as a Lyapunov function to satellite attitude tracking and robot navigation [37]. He showed that the dynamics arising from the natural motion of appropriately compensated mechanical systems are capable of integrating out the limit of a gradient system as well.
95 96 97 98 99 100 101 102 103 104 105 106 107 109
(38)
with the properties M T = M > 0, G T = −G, K T = K and A T = −A. The inertial matrix M is symmetric and positive definite, the gyric matrix G and the constraint damping forces A are skew-symmetric and the stiffness matrix K is symmetric. The generalized coordi˜ In state space form, Eq. (38) nates are represented by the vector q. becomes
x˜ =
a61 = −3 (c 2 c 5 + c 3 c 6 ) ,
60 62
81
80
Mq˜ + Gq˜ + Kq˜ = 0
2
58
61
The model can also be expressed in state-space form.
111 112 113
Using Eq. (33) in the cross-track part of Eq. (14) gives the following harmonically linearized cross-track equation
56
59
79
Consider a linearized equation of free motion with the general form [38,32]
3 2 c 5 + c 6 2 z˜ 4
˜
78
110
τ + 3c6 c5 sin τ cos τ + 3c5 c6 cos τ sin τ
+c 6 3 cos3 τ ≈
51
54
3
77
z˜ − a61 x˜ − a62 y˜ + (1 − a63 ) z˜ = b6
2
≈ (c 2 c 4 c 6 − c 3 c 4 c 5 ) + (c 2 c 6 − c 3 c 5 ) y˜ + c 2 2 + c 3 2 + c 4 2 z˜
46
53
(36)
5.1. Review of stability analysis
+c 5 c 4 sin τ + c 6 c 4 cos τ
45
50
y˜ + 2x˜ − a51 x˜ − a52 y˜ − a53 z˜ = b5
108
+4c 5 c 3 2 sin3 τ + 4c 3 2 c 6 sin2 τ cos τ
43
49
2
−8c 2 c 3 c 5 sin τ cos τ − 8c 2 c 3 c 6 sin τ cos τ
42
74
91
+c 3 c 5 sin τ cos τ + c 3 c 6 cos τ 1 2 ≈ c 2 + c 3 2 z˜ + 2 (c 2 c 5 + c 3 c 6 ) x˜
34
73
76
5. Application of total energy as a Lyapunov function for investigating stability of harmonic balance models
+c 5 c 2 2 sin3 τ + c 6 c 2 2 sin2 τ cos τ 2
72
75
(37)
z˜ x˜ = c 2 c 5 sin2 τ + c 2 c 6 cos τ sin τ + c 3 c 5 cos τ sin τ 1 +c 3 c 6 cos2 τ ≈ (c 2 c 5 + c 3 c 6 ) 2 z˜ x˜ 2 = 2c 2 c 3 c 5 sin2 τ cos τ + 2c 2 c 3 c 6 sin τ cos2 τ
33
37
71
The nonlinear terms in Eq. (32) are linearized as follows.
32
36
3
68 70
x˜ − 2 y˜ − (3 + a41 ) x˜ − a42 y˜ − a43 z˜ = b4
(c 2 c 6 − c 3 c 5 ) ,
2 3
67 69
The developed harmonic-balance model is summarized as follows
a51 = 3c 4 ,
11 12
5
0 I x˜ , −M−1 K −M−1 G
x˜ =
q˜ q˜
v˜ q˜ , q˜ =
1 T T 1 q˜ M q˜ + q˜ T Kq˜ 2 2
as a Lyapunov function.
115 116 117 118 119 120 121 122 123
(39)
124 125
A conservative gyric system in Eq. (39) is stable if it is statically stable. Stability can be inferred using the total system energy
(35)
114
126 127 128 129
(40)
130 131 132
JID:AESCTE
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.6 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
6
1
5.2. Investigation of the stability of harmonic balance model
is used as a Lyapunov function for the investigation of the stability ˜ x2 = q˜ the Lyapunov of harmonic balance model. Using x1 = q, function becomes
2
The harmonic balance model
3 4 5 6 7 8 9
x˜ = A(c)˜x + B(c)
(41)
is an inhomogeneous equation. Setting x˜ = 0 we have the equilibrium point given by the vector −1
10
x¯ = −[A(c)]
11
˜ then Using small perturbation parameter m,
12 13 14 15 16 17 18 19 20 21 22 23
B(c)
(42)
˜ = −[A(c)]−1 B(c) + m ˜ x˜ = x¯ + m
(43)
Upon substitution of Eq. (43) into Eq. (41) we have the homogeneous system
˜ ˜ = A(c)m m
(44)
˜ = 0. Clearly, the which has its equilibrium point at the origin m nature of the equilibrium point for Eq. (44) and Eq. (41) are iden˜ ˜ = A(c)m tical. Thus, we can apply Lyapunov stability analysis to m in the same way as for x˜ = A(c)˜x. The dynamics of the harmonic balance model may be written as
v˜ q˜ , q˜ =
1 2
x1T
x2T
KHB 0
0 1
x1 x2
(54)
This simplifies to
v˜ q˜ , q˜ = 0
(56)
Therefore, since the Lyapunov candidate function, v˜ q˜ , q˜ , is locally positive definite and its time derivative is zero, then the ˜ = 0, is proven to be Lyapunov stable. equilibrium point, m
k31 = −a61 , k32 = −a62 , k33 = − (−1 + a63 ) In Eq. (46), K H B is symmetric and G is skew symmetric. The Jacobian matrix of the conservative gyric system is
x˜ 0 = Ab(x0 )
28 29 30 31 32
and as a state space equation we have
d
dt
q˜ q˜
where,
33 34 35
KHB
= ⎡
q˜
(46)
− K H B q˜ − G q˜ − B
k11
k12
k13
⎤
⎡
0
g
0
⎤
= ⎣ k21 k22 k23 ⎦ , G = ⎣ − g 0 0 ⎦ , g = 2 k31
36 37
(45)
k32
0
k33
0
(47)
0
and
38 39 40
43 44 45 46 47
k11 = − (3 + a41 ) , k12 = −a42 , k13 = −a43
A=
48 49 50 51
0
1
−K H B
−G
(48)
(49)
and the linearized equation becomes
d
q˜ q˜
=
0 −K H B
1 −G
q˜ q˜
52
dt
53
which has a characteristic equation
(50)
56 57 58 59
62 63 64 65 66
77 78 79 80 82 83 84 85 86
(57)
110
Linear propagation of the calibrated initial condition using the traditional linear model has been seen to be more accurate than linear propagation of the true initial condition. For the purposes of the harmonic-balance model, this calibrated solution can be used as the generating solution. In forming the harmonic-balance model, this involves evaluating c(˜x0 ) as a function of the calibrated initial condition. The propagation of the harmonic-balance model, however, still uses the true initial condition, x0 .
112
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 111 113 114 115 116 117 118 119 120
λ2 + G λ + K H B = 0
(51)
and solutions
λ=
−G ±
G2
− 4K H B
2
60 61
76
90
54 55
75
89
42
27
q˜ + G q˜ + K H B q˜ − B = 0
74
88
k21 = −a51 , k22 = −a52 , k23 = −a53
26
73
87
6. Numerical simulations
41
25
71
81
The harmonic-balance method is essentially a two-step linearization approach. The nonlinear system is first linearized about the state-space origin, and the solution to the resulting linear system is used to define a generating solution. The nonlinear system is then linearized about the generating solution. The motivation is that the assumption of close proximity to the generating solution may be more accurate than the original assumption of close proximity to the origin. Recently, a method of coordinate calibration has been developed to extract more accurate linearized solutions for nonlinear systems [39]. The method takes advantage of coordinate transformations with alternative coordinates that describe the same system but enjoy a lower level of nonlinearity. Here, the calibration of the Cartesian coordinates, x, for the relative motion will be considered using transformations with the orbital-element differences, δ e. The nonlinear coordinate transformation, δ e = b(x), from Cartesian coordinates to orbital-element differences and the linearized coordinate transformation, x = Aδ e, from orbital element differences to Cartesian coordinates are used [40]. A calibrated initial condition is calculated as follows.
24
69
72
1 T 1 T 1 T 1 v q˜ , q = q˜ 1q˜ + q˜ 1q˜ + q˜ K H B q˜ + q˜ T K H B q˜ 2 2 2 2 T 1 T 1 = −G q˜ − K H B q˜ 1q˜ + q˜ 1 −G q˜ − K H B q˜ (55) 2 2
˜
68 70
>0
which is a radially positive definite function. Taking the unbounded derivative of v q˜ , q˜ gives
˜
67
(52)
Due to the nature of eigenvalues then the system is stable. Therefore, the system is globally exponentially stable because the poles of the system will always have negative real parts. The total energy of the system given as
1 T 1 T v˜ q˜ , q˜ = q˜ 1 q˜ + q˜ K H B q˜ 2 2
(53)
Case 1: The equations of motion in radial, along-track and cross-track directions are integrated using ode45. The orbital elements of the chief and deputy are shown in Table 1 and the orbit propagation was performed for six orbits. The simulations are carried out using two different initial conditions; true and CW calibrated initial conditions.
121 122 123 124 125 126 127 128
(a) Uncalibrated Generating Solution In this simulation, the CW, HB, cubic and nonlinear models are propagated using true initial conditions, and the vector of constants c is also evaluated using true initial conditions. The simulation results are illustrated in
129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.7 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
7
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52 54 55 56
Table 1 Chief and Deputy orbital elements for a spacecraft in LEO orbit.
57
Orbital elements
Chief
Deputy
58
Semimajor axis, a (km) Eccentricity, e Inclination, i (deg) RAAN, (deg) Argument of Perigee, ω (deg) True Anomaly, f (deg)
7500 0 4 5 10 25
7500 0.0012 4 5.01 15 20
59 60 61 62 63 64 65 66
118
Fig. 3. Propagated solutions for Example 1 using uncalibrated generating solution.
53
Fig. 3. These results indicate that the HB model approximates the nonlinear dynamics more accurately than the CW model, at least
119
over some initial time interval. To emphasize this another plot of approximation errors, focusing on the first 1000 seconds is shown. To quantify the approximation errors, the following error metric is defined for the average error.
ε(t ) = ρˆ (t ) − ρ (t ), ε¯ =
1 t sim
t sim
ε(t )dt
120 121 122 123 124 125
(58)
126 127
0
128
where ρ (t ) is the nonlinear solution for the relative position, ρˆ (t ) is one of the approximate linear solutions, and t sim is the length of time that has been simulated. The average error for each linear solution is shown in Table 2.
129 130 131 132
JID:AESCTE
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.8 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
8
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114 115
49
Fig. 4. Propagated solutions for Example 1 using calibrated generating solution.
50
116 117
51 52 53 54
Table 2 Errors over 1000 seconds for Example 1 using uncalibrated generating solution.
Table 3 Errors over 1000 seconds for Example 1 using calibrated generating solution.
118 119 120
55
Model
Average error (km)
Model
Average error (km)
121
56
CW HB
0.325 0.040
CW HB
0.962 0.0484
122
57
125
59 60 61 62 63 64 65 66
123 124
58
(b) Calibrated Generating Solution In this simulation, the HB, cubic and nonlinear models are propagated using true initial conditions, while the vector of constants c is evaluated using calibrated initial conditions. For purposes of a fair or complete comparison, the CW solution is propagated using the calibrated initial condition. The simulation results are illustrated in Fig. 4. The values for ε¯ for each linear model
is shown in Table 3. These results illustrate the improvement in accuracy of the HB model when using the calibrated generating solution, compared with using the uncalibrated generating solution. When using the uncalibrated generating solution, the HB model exhibited erroneous drift. When using the calibrated generating solution, the HB model maintained better accuracy over a larger time span. However, averaged over the simulation interval, the HB
126 127 128 129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.9 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
9
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116 117
51 52
Fig. 5. Propagated solutions for Example 2 using uncalibrated generating solution.
54 55 56 57 58 59
Table 4 Errors over 1000 seconds for Example 2 using uncalibrated generating solution. Model CW HB
0.325 0.040
61 63 64 65 66
ing the calibrated generating solution is able to provide good long term error without having to introduce initial error, unlike the calibrated CW solution.
model using calibrated generating solution was less accurate than the calibrated CW solution. The tradeoff, though, is the initial error introduced in the calibrated CW solution. The HB model has no such initial error. This is illustrated in the initial approximation errors. The HB model us-
120 121 122 123
Average error (km)
60 62
118 119
53
Case 2: A second example is constructed with initial conditions defined identical to Table 1 except the deputy’s semi-major axis is changed to 7505 km. The CW, HB, cubic and nonlinear models are propagated using true initial conditions, and the vector of constants c is also evaluated using true initial conditions. Fig. 5 shows the propagated solutions using the uncalibrated generating solution. The HB method provides good accuracy initially, but the error appears to grow rapidly later in the simulation
124 125 126 127 128 129 130 131 132
JID:AESCTE
10
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.10 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51 52 53 54 55 56 57
Table 5 Errors over 1000 seconds for Example 2 using calibrated generating solution.
59 61 62 63 64 65 66
(a) Calibrated generating solution The HB, cubic and nonlinear models are propagated using true initial conditions, while the vector of constants c is evaluated using calibrated initial conditions. The CW solution is propagated using the calibrated initial condition. Fig. 6 illustrates drifting orbit propagated solutions using calibrated generating solution and Table 5 shows errors over 1000 seconds.
118 119
interval. This may be related to how the HB model tries to approximate the nonlinear system as a linear, time varying system. To focus on the initial accuracy, the errors over the first 1000 seconds are shown in Table 4.
58 60
117
Fig. 6. Propagated solutions for Example 2 using calibrated generating solution.
Model
Average error (km)
CW HB
0.964 0.048
120 121 122 123 124 125 126
7. Conclusion
127 128
The harmonic-balance model of relative motion between chief spacecraft, in circular orbit, and deputy spacecraft has been presented. It results in an initial-condition dependent, linear model for the relative motion. The model is time-invariant. This paper has
129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.11 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
1 2 3 4 5 6 7 8
shown that harmonic-balance model of spacecraft relative motion gives a more accurate linearized model because it captured the motion better than the conventional method, Clohessy-Wiltshire model. The results obtained from the numerical simulations show that harmonic-balance model has lesser error than the ClohessyWiltshire model of satellite relative motion. Also, from the stability analysis using total energy as a Lyapunov function, we obtained that the harmonic balance model is Lyapunov stable.
9 10
8. Acknowledgment
11 12 13 14
The author would like to acknowledge Dr. Andrew Sinclair of Air Force Research Lab and Dr. Subhash Sinha of Auburn University for explaining the concept.
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] G. Hill, Researches in lunar theory, Am. J. Math. 1 (1878) 5–26. [2] W. Clohessy, R. Wiltshire, Terminal guidance system for satellite rendezvous, J. Guid. Control Dyn. 27 (9) (1960) 653–658. [3] S. Lee, I. Hwang, Reachable set computation for spacecraft relative motion with energy-limited low-thrust, Aerosp. Sci. Technol. 77 (2018) 180–188. [4] V. Pescea, R. Opromollab, S. Sarnoc, M. Lavagnaa, M. Grassib, Autonomous relative navigation around uncooperative spacecraft based on a single camera, Aerosp. Sci. Technol. 84 (2019) 1070–1080. [5] R. Ambrose, I. Nesnas, F. Chandler, B. Allen, T. Fong, L. Matthies, R. Mueller, Nasa technology roadmaps: ta 4: Robotics and autonomous systems, Technical report, NASA. [6] Y. Wang, H. Ji, Integrated relative position and attitude control for spacecraft rendezvous with ISS and finite-time convergence, Aerosp. Sci. Technol. 85 (2019) 234–245. [7] J. Wertz, Orbit and Constellation Design and Management, Microcosm Press and Springer, New York, 2009. [8] M. Macdonald, V. Badescu, The International Handbook of Space Technology, Springer-Verlag, Heidelberg, New York, 2014. [9] H. Curtis, Orbital Mechanics for Engineering Students, 2nd edition, Elsevier Ltd, 2010. [10] H. Schaub, J. Junkins, Analytical Mechanics of Space Systems, 3rd edition, AIAA Education Series, 2014. [11] G. Gomez, M. Marcote, Lindstedt-Poincaré solutions of Clohessy-Wiltshire equations, in: AAS/AIAA Astrodynamics Specialists Conference, 2005, pp. 1–13, Paper AAS 05-359. [12] D. Lee, G. Vukovich, Kinematically coupled spacecraft relative motion without attitude synchronization assumption, Aerosp. Sci. Technol. 45 (2015) 316–323. [13] R. Battin, Astronautical Guidance, McGraw-Hill, New York, 1964. [14] R. Battin, An Introduction to the Mathematics and Methods of Astrodynamics, AIAA Education Series, New York, 1987.
11
[15] T. Carter, M. Humi, Fuel-optimal rendezvous near a point in general Keplerian orbit, J. Guid. Control Dyn. 10 (6) (1987) 567–573. [16] H. London, Second approximation to the solution of the rendezvous equations, AIAA J. 1 (1963) 1691–1693. [17] Y. Zhang, L. Yang, Y. Zhu, H. Huang, W. Cai, Modeling and analysis of dynamics for spacecraft relative motion actuated by inter-satellite non-contacting force, Aerosp. Sci. Technol. 43 (2015) 236–244. [18] J. Li, C. Gao, C. Li, W. Jing, A survey on moving mass control technology, Aerosp. Sci. Technol. 82–83 (2018) 594–606. [19] J. Sullivan, S. Grimberg, S. DAmico, Comprehensive survey and assessment of spacecraft relative motion dynamics models, J. Guid. Control Dyn. 40 (2017) 1837–1859. [20] I. Huntley, R. Johnson, Linear and Nonlinear Differential Equations, Ellis Horwood Series, 1983. [21] D. Jordan, P. Smith, Nonlinear Ordinary Differential Equations: An Introduction for Scientists and Engineers, Oxford University Press, 2011. [22] A. Nafeh, D. Mook, Nonlinear Oscillations, John Wiley and Sons, New York, 1979. [23] S. Yamgoue, On the harmonic balance with linearization for asymmetric single degree of freedom non-linear oscillators, J. Nonlinear Dyn. 69 (2012) 1051–1062. [24] A. Belendez, E. Gimeno, M.L. Alvarez, M. Yebra, D. Mendez, Analytical approximate solutions for conservative nonlinear oscillators by modified rational harmonic balance method, Int. J. Comput. Math. (2020). [25] Z. Guo, A. Leung, H. Yang, Iterative homotopy harmonic balancing approach for conservative oscillator with strong odd-nonlinearity, J. Appl. Math. Model. 35 (2011) 1717–1728. [26] M. Ghadimi, H. Kaliji, Application of the harmonic balance method on nonlinear equations, World Appl. Sci. 22 (4) (2013) 532–537. [27] E. Grebenikov, Encyclopedia of Mathematics, Kluwer Acad. Publrs, ISBN 1402006098, 2011. [28] A. Ogundele, Nonlinear dynamics and control of spacecraft relative motion, Ph.D. thesis, Auburn University, Auburn, Alabama, 2017. [29] A. Ogundele, A. Sinclair, S. Sinha, Developing harmonic balance model of spacecraft relative motion, in: AAS/AIAA Spaceflight Mechanics Meeting, Napa, CA. [30] K. Huseyin, Vibrations and Stability of Multiple Parameter Systems, Sijthoff and Noordhoff International Publishers, The Netherlands, 1978. [31] W. Thompson, G. Tait, A Treatise on Natural Philosophy, 1st edition, Cambridge University Press, Cambridge, England, 1879. [32] W. Haddad, V. Chellaboina, Nonlinear Dynamical System and Control: A Lyapunov Based Approach, Princeton University Press, 2008. [33] N. Butenin, Elements of the Theory of Nonlinear Oscillations, Blaidsdell Publishing Company, New York, 1965. [34] R. Shieh, Energy and variational principles for generalized (gyroscopic) conservative problems, Int. J. Non-Linear Mech. 55 (1971) 495–509. [35] K. Huseyin, R. Plaut, Transverse vibrations and stability of systems with gyroscopic forces, J. Struct. Mech. 3 (2) (1974) 163–177. [36] K. Huseyin, R. Plaut, Standard forms of the eigenvalue problems associated with gyroscopic systems, J. Sound Vib. 45 (1) (1976) 29–37. [37] D. Koditschek, The application of total energy as a Lyapunov function for mechanical control systems, Contemp. Math. 97 (1989) 131–157. [38] P. Hughes, Spacecraft Attitude Dynamics, Dover Publication, Inc., New York, 2004. [39] A. Sinclair, R. Sherrill, T. Lovell, Calibration of linearized solutions for spacecraft relative motion, J. Guid. Control Dyn. 37 (2014) 1362–1367. [40] K. Alfriend, H. Schaub, D. Gim, Gravitational perturbations, nonlinearity and circular orbit assumption effects on formation flying control strategies, in: AAS Rocky Mountain Guidance and Control Conference, vol. 104, 2000, pp. 139–158.
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65
131
66
132
[m5G; v1.261; Prn:5/02/2020; 15:29] P.1 (1-11)
Aerospace Science and Technology ••• (••••) ••••••
1
Contents lists available at ScienceDirect
67 68
2 3
Aerospace Science and Technology
4
69 70 71
5
72
6
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13
Modeling and analysis of nonlinear spacecraft relative motion via harmonic balance and Lyapunov function ✩
16 17
81
Mechanical and Aerospace Engineering Department, Naval Postgraduate School, Monterey, CA, United States of America
83
82 84 85
19
a r t i c l e
i n f o
a b s t r a c t
86 87
21 22 23 24 25 26 27 28 29 30 31 32
79
Ayansola D. Ogundele
18 20
78 80
14 15
77
Article history: Received 14 October 2019 Received in revised form 27 October 2019 Accepted 1 February 2020 Available online xxxx Communicated by Efstathios Bakolas Keywords: Relative motion Cubic approximation Harmonic balance Lyapunov function Total energy
The harmonic balance method, an approximate method for the study of nonlinear oscillating systems described by ordinary nonlinear differential equations, is applied to the cubic polynomial approximation of the relative motion to develop approximate periodic solutions of the motion. Firstly, the cubic nonlinear terms neglected in the formulation of linear time invariant Clohessy-Wilshire (CW) equations are obtained from the original nonlinear equation of motion. Secondly, using the CW solutions as the generating solutions, harmonic balance technique is applied to the cubic nonlinear terms and correction terms are obtained. The new and improved harmonically linearized model contains CW equations and the correction terms. Using total energy as a Lyapunov function, the stability of the new model is determined. The additional terms in the new harmonically linearized spacecraft relative motion enabled us to obtain excellent approximate periodic solutions of the nonlinear equations as evidenced from the numerical simulation results. © 2020 Elsevier Masson SAS. All rights reserved.
88 89 90 91 92 93 94 95 96 97 98
33
99
34
100
35
101
36
102
37 38
1. Introduction
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
George Hill (1878) in the nineteenth century discovered that a great insight into the dynamics of nearby objects can be gained by linearizing the relative equations of motion and ignoring the longterm effects of the natural perturbations in the context of celestial mechanics [1]. In 1960, W. R. Clohessy and R. S. Wiltshire used this approach when researching on terminal guidance for satellite rendezvous in the context of astrodynamics [2]. The relative motion dynamics between two spacecraft, for example Chief and Deputy spacecraft, are generally analyzed using Clohessy-Wilshire (HCW) equations. In this model, a reference frame centered on a reference spacecraft, called Chief spacecraft, that orbits the Earth in a circular trajectory is defined. The assessment of the operational capability of the spacecraft in formation flying can be done using the knowledge on a reachable set for their relative motion [3]. Proper modeling of relative motion dynamics is desirable as this will enhance spacecraft autonomy which is required by the next generation of spacecraft [4,5]. It was established that [6] the relative position dynamics is more suitable for actual measurement circumstances and applicable for arbitrary orbital forms.
59 60 61 62 63 64 65 66
✩
This document is a revised version of the paper presented at 2016 AAS/AIAA Spaceflight Mechanics Meeting. E-mail address: [email protected]. https://doi.org/10.1016/j.ast.2020.105761 1270-9638/© 2020 Elsevier Masson SAS. All rights reserved.
The motion of Deputy spacecraft moving relative to the Chief is studied from a reference frame fixed at the center of the reference spacecraft called Local Vertical Local Horizontal frame [2,7–12]. The orbital radius of the reference spacecraft is assumed to be much larger than the relative separation between the spacecraft. Over a long time this approach has significant errors. Because of this limitation there have being some modifications in one way or the other over a period of time. Battin uses universal variables which allows it to be used for elliptical and hyperbolic orbits [13,14]. Also, Carter and Humi modified the dynamics for eccentric reference orbits [15]. London [16] developed and solved second order polynomial approximation of satellite equations of motion but a better accuracy can be obtained using higher order polynomial. In [17] detailed modeling and analysis of spacecraft relative motion dynamics actuated by inter-satellite non-contacting force is given. Survey of the technological difficulties and developmental perspectives in spacecraft dynamic analysis are highlighted in [18]. A comprehensive survey and assessment of spacecraft relative motion can be found in [19]. The harmonic balance (HB) method, or harmonic linearization, is a powerful technique for the analysis of high-frequency nonlinear oscillators and for obtaining analytic solutions to nonlinear ordinary differential equations describing oscillatory systems [20,21]. The method became prominent in 1990s because of its ability to calculate and capture the steady-state behavior of nonlinear systems directly and high applicability. Harmonic balance
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.2 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
method is based on trigonometric functions and it is possible to formulate this procedure using any complete set of periodic functions such as Jacobi elliptic functions. It can be adapted to construct a linear substitute for the original equation. The difficulties in using this method are that calculating the amplitudes and the angular frequency may become algebraically intensive and existing formulations of the harmonic balance procedure do not allow it to be applied to non-conservative oscillators because these types of solutions have solutions involving transient behaviors. In order to obtain a consistent solution using the method one needs to either know a great deal about the solution a priori or to carry enough terms in the solution and check the order of the coefficients of all the neglected harmonics, otherwise one might obtain inaccurate approximation [22,23]. Before applying harmonic balance the governing equations are firstly linearized to a lower order polynomial so as to prevent solving complicated nonlinear algebraic equations. This leads to linear algebraic equations instead of nonlinear one. The harmonic linearization has been applied successfully to several nonlinear systems [24–26]. The main contributions of this paper are itemized as follows:
22 23
• A new harmonically linearized spacecraft relative motion is
24 25
•
26 27 28
•
29 30 31
•
32 33 34
•
developed making use of assumed generating solution. The method is amenable based on the fact that it can be applied to solve differential equations whose nonlinear terms are not small. The resulted harmonic balance linearized form gave a better approximation of periodic solutions of the original nonlinear term. The new model is useful for analyzing spacecraft formation flying, rendezvous and proximity operations. Using total energy as a Lyapunov function the developed harmonic balance model stability was investigated.
35 36 37 38 39 40 41 42 43 44 45 46 47 48
In this paper, we developed harmonically linearized spacecraft relative motion using harmonic balance technique. To the best of our knowledge, this method is rarely used in the papers dealing with the spacecraft relative motion dynamics. The approach is a two stage procedure. The first approach involves linearization of the spacecraft nonlinear equation of motion into a cubic polynomial model using Binomial theorem and the second approach involves the application of harmonic linearization technique into the cubic model. This leads to the formulation of new harmonic balance relative equation of motion consisting of improved, three linear, second order ODEs. We investigated the stability of the model using total energy as Lyapunov function and it was found to be stable.
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
2. Review of harmonic balance method The review of harmonic balance method is adapted from the References [27,28,22,20,21,29]. This method is most commonly used to study the nonlinear vibration problems. It may be used to determine the approximate periodic and quasi-periodic oscillations, periodic and quasi-periodic conditions in automatic control theory, as well as stationary conditions, and in the studies of their stability. If a periodic solution does exist, it may be sought in the form of a Fourier series, whose coefficients are determined by requiring the series to satisfy the equation of motion. Consider a nonlinear perturbed oscillator differential equation
x¨ + ω x + ε f (x, x˙ ) = 0, 2
(0 < ε 1)
(1)
where ε is a small parameter. Harmonic linearization technique involves the replacement of the nonlinear forcing function F (x, x˙ ) = ε f (x, x˙ ) by the equivalent linear force [27]
67 68 69 70 71 72 73 74 75 76 77 78 79
Fig. 1. Steps followed in finding an approximate model for the nonlinear equation using Harmonic Balance method.
80 81 82
F l (x, x˙ ) = kx + λ˙x
(2)
83 84
The parameters k and λ can be computed from
ε ∫ f (a cos ψ, −aω sin ψ) cos ψ dψ πa 0 ε 2π ∫ f (a cos sψ, −aω sin ψ) sin ψ dψ λ (a) = − π aω 0
k (a) =
85 86
2π
87 88
(3)
where ψ = ωt + θ . If x = acos (ωt + θ), a= constant, ω = constant, the nonlinear force F (x, x˙ ) is a periodic function of time, and its Fourier series expansion contains an infinite number of harmonics, having the frequencies nω , n = 1, 2, . . .. The forcing function F (x, x˙ ) has the form
F (x, x˙ ) =
∞
89 90 91 92 93 94 95 96 97
F n cos (nωt + θn )
(4)
n =0
98 99 100
The term F 1 cos(ωt + θ1 ) is called the fundamental harmonic of the expansion. The amplitude and the phase of the linear function F l coincide with the respective characteristics of the fundamental harmonic of the nonlinear force. For the differential equation
101 102 103 104 105
x¨ + ω2 x + F (x, x˙ ) = 0
(5)
106 107
which is typical in the theory of quasi-linear oscillations, the harmonic balance method consists in replacing F (x, x˙ ) by the equivalent linear force to give
x¨ + λ˙x + k1 x = 0
108 109 110 111
(6)
112
where k1 = k + ω2 , λ is the equivalent damping coefficient and k1 is the equivalent elasticity coefficient. In the harmonic balance method the frequency of the oscillation depends on the amplitude (through the quantities k and λ). By traditional linearization about the origin we have
113
x¨ + ω x = 0
119
2
(7)
Comparing the harmonic linearization result in Eq. (6) with Eq. (7) it can be easily seen that harmonic linearized equation has correction terms. In this work, harmonic balance method is adapted to construct a new and improved linearized model of the spacecraft relative motion. Fig. 1 illustrates the steps followed in finding an approximate model for a nonlinear system. In developing the harmonicbalance model for spacecraft relative motion, polynomial approximation is applied to the nonlinear equations in LVLH components and it produced cubic equations of motion containing linear and nonlinear terms. The linear equations correspond to the ClohessyWiltshire equations. For the harmonic linearization, the nonlinear
114 115 116 117 118 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.3 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
3
The LHS of Eq. (10) shows the differentiation in inertial frame while that in the RHS shows the differentiation in Hill frame. Therefore, differentiating Eq. (8) twice and using Eq. (10) gives the kinematics equation of motion as
1 2 3 4
¨ + ρ¨ + 2ω × ρ˙ + ω × (ω × ρ ) r¨ = R
6 7 8
(11)
For the case with the chief in circular orbit the mean motion is ˙f = n = μ/ R 3 . This reduces Eq. (11) to
9 10
11
r¨ = x¨ − 2ny − n2 x −
12
μ ˆ
i + y¨ + 2nx˙ − n2 y ˆj + z¨ kˆ
R2
(12)
Using Eq. (12) in the second part of Eq. (9) gives the second order nonlinear differential equations of relative motion as
14 15 16
x¨ − 2ny − n2 x −
17 18
μ R2
y¨ + 2nx˙ − n y = −
20
Fig. 2. Relative Motion of Deputy spacecraft with respect to Chief Spacecraft.
z¨ = −
22
25 26
terms are evaluated along the generating solution to produce linear correction terms. The correction terms are combined with the original linear equations (CW equations) to produce the harmonic balance model of the spacecraft relative motion.
27 28
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
˜
58 59 60 61 62 63 64 65 66
72 73 74 75 76 77 78
horizontal frame (xyz) with unit vectors
ˆi, ˆj, kˆ , referred to as
Hill frame, is attached to the chief and it is rotating with it. The x axis is aligned along the outward radial direction, the z axis is aligned with the direction perpendicular to the orbit plane and it is parallel to the orbit momentum vector, and the y axis is normal to the other two axes and completes the right-hand coordinate system. The in-plane motion is defined by x and y while the out-of-plane motion is defined by the z axis. Expressed in the Hill frame, the deputy’s position relative to the chief is denoted as ρ = xˆi + yˆj + zkˆ = (x, y , z)T . The unperturbed rotating Hill (LVLH) frame relative to the inertial frame has constant angular velocity ˙ = ¨f kˆ = 0. From of ω = ˙f kˆ = Rh2 kˆ and angular acceleration of ω Fig. 1, R is the chief spacecraft position vector relative to the Earth and can be expressed as R = R ˆi in the Hill frame. The angular mo˙ = R 2 ˙f kˆ and its mentum expressed in the Hill frame is h = R × R 2˙ magnitude is h = R f . The deputy spacecraft position vector r relative to the center of the gravitational field is
r=R+ρ
r3
˜
(8)
˜
z + z˜ = V x˜ , y˜ , z˜
(13)
88
¨ =− R r¨ = −
R3
(14)
μ r3
r
(9)
2 Nd ρ
dt 2
= ρ¨ + 2ω × ρ˙ + ω × (ω × ρ )
(10)
94 95 97
99
2
3 2
y˜ + 2
3
100
z˜ − 6x˜ y˜ − 6x˜ z˜ − 4x˜ 2
2
2 3 2 3 3 2 U x˜ , y˜ , z˜ = 3 y˜ x˜ − 6 y˜ x˜ + y˜ z˜ − y˜ 2 2 3 3 V x˜ , y˜ , z˜ = 3 z˜ x˜ − 6 z˜ x˜ 2 + z˜ y˜ 2 − z˜ 3 2 2
2
3
101 102
(15)
105 106 107 108 109 110 111 112
x˜ − 2 y˜ − 3x˜ = 0 y˜ + 2x˜ = 0
103 104
Eq. (14) is the nonlinear cubic approximation model of spacecraft relative motion with chief in circular orbit. For ε = 0, that is neglecting the nonlinear terms, the system in Eqs. (14) reduced to the linear equations
113
(16)
114 115
˜
z + z˜ = 0
116
These equations are the well-known normalized Clohessy-Wilshire (CW) Equations. They are coupled equations of motion in x and y. Using linear theory, the homogeneous equations have the solutions [9,29]
118
117 119 120 121 122 123 124
(17)
127 128
τ = nt and
c 1 = 2x˜ 0 + y˜ 0 , c 2 = x˜ 0 , c 3 = −3x˜ 0 − 2 y˜ c 4 = −2x˜ 0 + y˜ 0 , c 5 = z˜ 0 , c 6 = z˜ 0
125 126
z˜ = c 5 sin τ + c 6 cos τ where,
Using transport theorem we obtain
91
96
y˜ = −3c 1 τ + 2c 2 cos τ − 2c 3 sin τ + c 4
R
90
93
x˜ = 2c 1 + c 2 sin τ + c 3 cos τ
μ
89
92
The Keplerian acceleration of the chief and deputy spacecraft are
86
98
T x˜ , y˜ , z˜ = −3x˜ +
85 87
z
where,
81 83
y
r3
μ
80
84
μ
y + 2x = U x˜ , y˜ , z˜
In this section, the cubic polynomial approximation model of deputy spacecraft motion relative to the chief spacecraft in circular orbit is derived. The derivation follows the approach in [10,29,28]. Fig. 2 shows relative motion of deputy spacecraft with respect to the chief spacecraft in circular orbit. A local-vertical/local
( R + x)
r3
x˜ − 2 y˜ − 3x˜ = T x˜ , y˜ , z˜
3. Nonlinear cubic approximation model of relative motion
56 57
70
82
μ
Applying Binomial theorem expansion to Eq. (13) and normalizing the resulting equations we have
29 30
=−
2
19
24
69
79
13
23
68
71
5
21
67
129 130
(18)
131 132
JID:AESCTE
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.4 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
4
1 2
4. Harmonic linearization of nonlinear cubic approximation model
Substituting the linear approximations in Eqs. (23) and (24) into Eq. (23) and use the results in the radial part of Eq. (14) gives
3 4 5 6 7 8 9
In this section, harmonic linearization is applied to the cubic approximation model of spacecraft relative motion using the bounded, steady-state solutions of CW equations as the generating solutions. In order to obtain bounded solutions, the secular term in the along-track is eliminated by setting c 1 = 0 and this leads to
10
x˜ = c 2 sin τ + c 3 cos τ
11
y˜ = 2c 2 cos τ − 2c 3 sin τ + c 4
12 13
z˜ = c 5 sin τ + c 6 cos τ
14 15 16 18 19 20 21 22 23
The radial nonlinear part, from Eq. (14), is given by
T x˜ , y˜ , z˜ = −3x˜ +
cos2 τ =
26 27 28
2
y˜ + 2
3
z˜ − 6x˜ y˜ − 6x˜ z˜ − 4x˜ 2
2
2
2
3
cos
3
τ=
1 2 1 4
1
sin2 τ =
(1 + cos 2τ ),
2 3
(3 cos τ + cos 3τ ),
sin
τ=
1
(3 sin τ − sin 3τ )
4
31 32
(21)
33 34 35
The nonlinear terms are linearized as follows
x˜ = 2c 2 c 3 sin τ cos τ 2
+ c 22 sin2
τ
+ c 32 cos2
τ
Substituting Eq. (21) into Eq. (22) and eliminating the higher order harmonic of time yields
41
Similar simplification is performed for the remaining terms as
38 39
42 43 44 45 46 47 48
y˜ = 2
τ
51 52 53 54 55
≈
56
59 60 61 62 63 64
x˜ z˜
2
c 22
τ
τ
τ
τ
τ
τ
τ
τ
+ c 32
+ c 42
= c 3 c 62 cos3 +c 3 c 52 cos
τ
τ
U x˜ , y˜ , z˜ = 3 y˜ x˜ − 6 y˜ x˜ 2 +
3 2
y˜ z˜ 2 −
3
τ sin τ
+2c 3 c 5 c 6 cos τ sin τ 1 ≈ (2c 2 c 5 + 2c 3 c 6 ) z˜ + c 52 + c 62 x˜
4 x˜ 3 = c 33 cos3 τ + c 23 sin3 τ + 3c 2 c 32 cos2 τ sin τ
3 2 ≈ c 2 + c 32 x˜ 4
80 81 82 83 84 86
2
88 89
(27)
y˜ x˜ = 4c 2 c 3 sin τ cos
2
2
3
2
95 97
100 101
3
τ + 2c2 sin τ cos τ + 2c2 c3 cos τ
102 103 104
−2c 3 sin τ cos τ + 2c 2 c 3 c 4 sin τ cos τ
105
+c 2 2 c 4 sin2 τ + c 3 2 c 4 cos2 τ 1 2 = c 2 c 4 + c 3 2 c 4 + c 2 2 + c 3 2 y˜ 2
94
99
2
3
93
98
The remaining nonlinear terms are harmonically linearized as 2
92
96
(28)
−2c 3 sin τ cos τ + c 2 c 4 sin τ + c 4 c 3 cos τ ≈ c 4 x˜ 2
90 91
2
106 107 108 109 110
3
+2c 2 c 6 cos τ − 2c 3 c 5 sin τ
111
−4c 3 c 5 c 6 sin2 τ cos τ − 2c 3 c 6 2 sin τ cos2 τ
113
2
+c 4 c 6 2 cos2 τ 1 2 = c 5 + c 6 2 c 4 + c 5 2 + c 6 2 y˜ 4
+4 (c 2 c 6 − c 3 c 5 ) z˜ y˜ 3 = 2c 2 c 4 2 cos τ + 8c 2 2 c 4 cos2 τ + 2c 2 c 4 2 cos τ 2
+8c 2 c 4 cos τ − 8c 2 c 3 c 4 sin τ cos τ
2
79
y˜ 3
+c 4 c 5 sin τ + 2c 4 c 5 c 6 sin τ cos τ
(24)
77 78
Harmonic linearization of Eq. (27) is performed as follows. Substituting the trigonometric identities in Eq. (21) into Eq. (27), eliminating the higher order harmonics and rearrange the result in the form of assumed linear solution yields
2
τ
τ sin τ + 2c2 c5 c6 cos τ sin2 τ
+3c 22 c 3 cos τ sin2 τ
+ 3c 62
76
87
2
+ c 2 c 62 cos2
+ 3c 52
The along-track nonlinear part, from Eq. (14), is
2
τ
x˜
+ c 2 c 52 sin3 2
+ 18c 42
4 y˜ z˜ 2 = 2c 2 c 5 2 cos τ sin2 τ + 4c 2 c 5 c 6 sin τ cos2 τ
τ
τ
65 66
(23)
2 + c 3 c 42 cos τ
+c 2 c 42 sin + 4c 22 c 3 cos3 + 4c 2 c 32 sin3 +4c 22 c 4 cos sin − 4c 32 c 4 cos sin −8c 2 c 32 cos2 sin − 8c 22 c 3 cos sin2 +4c 2 c 3 c 4 cos2 − 4c 2 c 3 c 4 sin2
+ 6c 32
3
τ − 4c3 c4 sin τ + 4c2 c4 cos τ
x˜ y˜ 2 = 4c 23 cos2 τ sin τ + 4c 33 cos τ sin2 τ
50
58
τ
+ 4c 32 sin2
4
6c 22
(26)
−4c 2 c 3 2 sin2 τ cos τ − 2c 2 2 c 3 sin3 τ
−8c 2 c 3 cos τ sin τ + c 42 ≈ −2c 4 y˜ + 2c 22 + 2c 32 + 3c 42 1 2 z˜ 2 = c 62 cos2 τ + c 52 sin2 τ + 2c 5 c 6 cos τ sin τ ≈ c 5 + c 62
49
57
4c 22 cos2
75
85
(22)
40
37
1
74
y˜ x˜ = 2c 2 2 sin τ cos τ + 2c 2 c 3 cos2 τ − 2c 2 c 3 sin2 τ
1 2 1 2 x˜ 2 = c 2 c 3 sin 2τ + c 3 − c 22 cos 2τ + c 2 + c 32 2 2 1 2 2 ≈ c + c3 2 2
36
73
2 a42 = 3c 4 ,
b4 =
70 72
4.2. Harmonic linearization of along-track equation of motion
(1 − cos 2τ ),
29 30
−18c 22 − 18c 32 − 12c 42 − 3c 52 − 3c 62 ,
(20)
Harmonic linearization of each of the nonlinear terms is performed by substituting the generating solution in Eq. (19) into the nonlinear terms and using the following trigonometric identities:
24 25
3
2
1
Eq. (25) is the radial harmonically linearized relative motion. The additional terms in x, y , z and the nonhomogeneous term b4 are the radial corrections to the equation of motion. They are approximations of the cubic powers that were retained in the expansion of the nonlinear equations of motion.
4.1. Harmonic linearization of radial equation of motion
17
(25)
71
where,
a43 = −3 (c 2 c 5 + c 3 c 6 ), (19)
68 69
x˜ − 2 y˜ − (3 + a41 ) x˜ − a42 y˜ − a43 z˜ = b4
a41 =
67
−16c2 2 c 3 sin τ cos2 τ + 8c 2 3 cos3 τ 2
2
2
+8c 2 c 3 sin τ cos τ − 2c 3 c 4 sin τ −8c 2 c 3 c 4 sin τ cos τ + 8c 3 2 c 4 sin2 τ 2
2
+16c2 c 3 sin τ cos τ = 3c 2 2 c 4 + 3c 3 2 c 4 − 2c 4 3 + 3c 2 2 + 3c 3 2 + 3c 4 2 y˜
112
(29)
114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.5 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
1 2 3 4 5 6 7 8 9 10
Substituting these linear approximations into Eq. (27) and using the results in the along-track part of Eq. (14) leads to the harmonically linearized along-track equation of motion
The additional terms in x, y , z and the nonhomogeneous term b3 (c ) are the corrections to the cross-track direction.
y˜ + 2x˜ − a51 x˜ − a52 y˜ − a53 z˜ = b5
4.4. Model summary
13 14 15 16 17 18 19
(30)
where,
3 a52 = −16c2 2 − 16c3 2 − 12c4 2 + c 5 2 + c 6 2 , a53 = b5 =
8 3
(31)
8
−16c4 c 2 2 − 16c4 c 3 2 + c 4 c 5 2 + c 4 c 6 2 + 8c 4 3
The additional terms in x, y , z and the nonhomogeneous term b5 (c ) are the corrections to the along-track direction. 4.3. Harmonic linearization of cross-track equation of motion
20 21
The cross-track nonlinear part, from Eq. (14), is
22 23 24 25 26 27 28 29 30 31
V x˜ , y˜ , z˜ = 3 z˜ x˜ − 6 z˜ x˜ 2 +
2
z˜ y˜ 2 −
3 2
z˜ 3
(32)
2
2
35
3
4 z˜ y˜ 2 = 4c 5 c 2 2 cos2 τ sin τ + 4c 6 c 2 2 cos3 τ
38
(33) 2
39
+4c 2 c 4 c 5 cos τ sin τ + 4c 2 c 4 c 6 cos τ
40
−4c 3 c 4 c 5 sin2 τ − 4c 3 c 4 c 6 sin τ cos τ
41
2
44
2
47 48
z˜ = c 5 sin 3
3
52
55
2
2
z − a61 x˜ − a62 y˜ + (1 − a63 ) z˜ = b6
57
where,
2
(34)
63 64 65 66
82
⎡
⎤ ⎡ 0 0 0 1 0 0 0 0 0 0 0 1 0⎥ ⎢ 0 ⎢ ⎢ ⎢ 0 0 0 0 0 1⎥ 0 ⎢ ⎥ ⎢ x˜ = ⎢ + 3 + a41 (c) a42 (c) a43 (c) 0 2 0 ⎥ ⎢ b 4 ( c) ⎣ ⎦ ⎣ a51 (c) a52 (c) a53 (c) b 5 ( c) −2 0 0 a61 (c) a62 (c) a63 (c) − 1 0 0 0 b 6 ( c) x˜ = A (c) x˜ + B (c)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦
a62 =
3
(c 2 c 6 − c 3 c 5 ) 2 1 a63 = 12c 4 2 − 9c 5 2 − 9c 6 2 , 8 1 12c 4 2 − 9c 5 2 − 9c 6 2 b6 = 8
83 84 85 86 87 88 89 90 92 93 94
In this section, total energy is applied as a Lyapunov function to investigate the stability of newly developed harmonic balance model of spacecraft relative motion. Several authors have studied vibrations and stability of gyroscopic conservative systems [30–32]. General discussions and developments on gyroscopic systems, and a variety of investigations ranging from important practical problems to new general theories can be found in [33–36]. In 1989, Koditschek applied total energy as a Lyapunov function to satellite attitude tracking and robot navigation [37]. He showed that the dynamics arising from the natural motion of appropriately compensated mechanical systems are capable of integrating out the limit of a gradient system as well.
95 96 97 98 99 100 101 102 103 104 105 106 107 109
(38)
with the properties M T = M > 0, G T = −G, K T = K and A T = −A. The inertial matrix M is symmetric and positive definite, the gyric matrix G and the constraint damping forces A are skew-symmetric and the stiffness matrix K is symmetric. The generalized coordi˜ In state space form, Eq. (38) nates are represented by the vector q. becomes
x˜ =
a61 = −3 (c 2 c 5 + c 3 c 6 ) ,
60 62
81
80
Mq˜ + Gq˜ + Kq˜ = 0
2
58
61
The model can also be expressed in state-space form.
111 112 113
Using Eq. (33) in the cross-track part of Eq. (14) gives the following harmonically linearized cross-track equation
56
59
79
Consider a linearized equation of free motion with the general form [38,32]
3 2 c 5 + c 6 2 z˜ 4
˜
78
110
τ + 3c6 c5 sin τ cos τ + 3c5 c6 cos τ sin τ
+c 6 3 cos3 τ ≈
51
54
3
77
z˜ − a61 x˜ − a62 y˜ + (1 − a63 ) z˜ = b6
2
≈ (c 2 c 4 c 6 − c 3 c 4 c 5 ) + (c 2 c 6 − c 3 c 5 ) y˜ + c 2 2 + c 3 2 + c 4 2 z˜
46
53
(36)
5.1. Review of stability analysis
+c 5 c 4 sin τ + c 6 c 4 cos τ
45
50
y˜ + 2x˜ − a51 x˜ − a52 y˜ − a53 z˜ = b5
108
+4c 5 c 3 2 sin3 τ + 4c 3 2 c 6 sin2 τ cos τ
43
49
2
−8c 2 c 3 c 5 sin τ cos τ − 8c 2 c 3 c 6 sin τ cos τ
42
74
91
+c 3 c 5 sin τ cos τ + c 3 c 6 cos τ 1 2 ≈ c 2 + c 3 2 z˜ + 2 (c 2 c 5 + c 3 c 6 ) x˜
34
73
76
5. Application of total energy as a Lyapunov function for investigating stability of harmonic balance models
+c 5 c 2 2 sin3 τ + c 6 c 2 2 sin2 τ cos τ 2
72
75
(37)
z˜ x˜ = c 2 c 5 sin2 τ + c 2 c 6 cos τ sin τ + c 3 c 5 cos τ sin τ 1 +c 3 c 6 cos2 τ ≈ (c 2 c 5 + c 3 c 6 ) 2 z˜ x˜ 2 = 2c 2 c 3 c 5 sin2 τ cos τ + 2c 2 c 3 c 6 sin τ cos2 τ
33
37
71
The nonlinear terms in Eq. (32) are linearized as follows.
32
36
3
68 70
x˜ − 2 y˜ − (3 + a41 ) x˜ − a42 y˜ − a43 z˜ = b4
(c 2 c 6 − c 3 c 5 ) ,
2 3
67 69
The developed harmonic-balance model is summarized as follows
a51 = 3c 4 ,
11 12
5
0 I x˜ , −M−1 K −M−1 G
x˜ =
q˜ q˜
v˜ q˜ , q˜ =
1 T T 1 q˜ M q˜ + q˜ T Kq˜ 2 2
as a Lyapunov function.
115 116 117 118 119 120 121 122 123
(39)
124 125
A conservative gyric system in Eq. (39) is stable if it is statically stable. Stability can be inferred using the total system energy
(35)
114
126 127 128 129
(40)
130 131 132
JID:AESCTE
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.6 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
6
1
5.2. Investigation of the stability of harmonic balance model
is used as a Lyapunov function for the investigation of the stability ˜ x2 = q˜ the Lyapunov of harmonic balance model. Using x1 = q, function becomes
2
The harmonic balance model
3 4 5 6 7 8 9
x˜ = A(c)˜x + B(c)
(41)
is an inhomogeneous equation. Setting x˜ = 0 we have the equilibrium point given by the vector −1
10
x¯ = −[A(c)]
11
˜ then Using small perturbation parameter m,
12 13 14 15 16 17 18 19 20 21 22 23
B(c)
(42)
˜ = −[A(c)]−1 B(c) + m ˜ x˜ = x¯ + m
(43)
Upon substitution of Eq. (43) into Eq. (41) we have the homogeneous system
˜ ˜ = A(c)m m
(44)
˜ = 0. Clearly, the which has its equilibrium point at the origin m nature of the equilibrium point for Eq. (44) and Eq. (41) are iden˜ ˜ = A(c)m tical. Thus, we can apply Lyapunov stability analysis to m in the same way as for x˜ = A(c)˜x. The dynamics of the harmonic balance model may be written as
v˜ q˜ , q˜ =
1 2
x1T
x2T
KHB 0
0 1
x1 x2
(54)
This simplifies to
v˜ q˜ , q˜ = 0
(56)
Therefore, since the Lyapunov candidate function, v˜ q˜ , q˜ , is locally positive definite and its time derivative is zero, then the ˜ = 0, is proven to be Lyapunov stable. equilibrium point, m
k31 = −a61 , k32 = −a62 , k33 = − (−1 + a63 ) In Eq. (46), K H B is symmetric and G is skew symmetric. The Jacobian matrix of the conservative gyric system is
x˜ 0 = Ab(x0 )
28 29 30 31 32
and as a state space equation we have
d
dt
q˜ q˜
where,
33 34 35
KHB
= ⎡
q˜
(46)
− K H B q˜ − G q˜ − B
k11
k12
k13
⎤
⎡
0
g
0
⎤
= ⎣ k21 k22 k23 ⎦ , G = ⎣ − g 0 0 ⎦ , g = 2 k31
36 37
(45)
k32
0
k33
0
(47)
0
and
38 39 40
43 44 45 46 47
k11 = − (3 + a41 ) , k12 = −a42 , k13 = −a43
A=
48 49 50 51
0
1
−K H B
−G
(48)
(49)
and the linearized equation becomes
d
q˜ q˜
=
0 −K H B
1 −G
q˜ q˜
52
dt
53
which has a characteristic equation
(50)
56 57 58 59
62 63 64 65 66
77 78 79 80 82 83 84 85 86
(57)
110
Linear propagation of the calibrated initial condition using the traditional linear model has been seen to be more accurate than linear propagation of the true initial condition. For the purposes of the harmonic-balance model, this calibrated solution can be used as the generating solution. In forming the harmonic-balance model, this involves evaluating c(˜x0 ) as a function of the calibrated initial condition. The propagation of the harmonic-balance model, however, still uses the true initial condition, x0 .
112
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 111 113 114 115 116 117 118 119 120
λ2 + G λ + K H B = 0
(51)
and solutions
λ=
−G ±
G2
− 4K H B
2
60 61
76
90
54 55
75
89
42
27
q˜ + G q˜ + K H B q˜ − B = 0
74
88
k21 = −a51 , k22 = −a52 , k23 = −a53
26
73
87
6. Numerical simulations
41
25
71
81
The harmonic-balance method is essentially a two-step linearization approach. The nonlinear system is first linearized about the state-space origin, and the solution to the resulting linear system is used to define a generating solution. The nonlinear system is then linearized about the generating solution. The motivation is that the assumption of close proximity to the generating solution may be more accurate than the original assumption of close proximity to the origin. Recently, a method of coordinate calibration has been developed to extract more accurate linearized solutions for nonlinear systems [39]. The method takes advantage of coordinate transformations with alternative coordinates that describe the same system but enjoy a lower level of nonlinearity. Here, the calibration of the Cartesian coordinates, x, for the relative motion will be considered using transformations with the orbital-element differences, δ e. The nonlinear coordinate transformation, δ e = b(x), from Cartesian coordinates to orbital-element differences and the linearized coordinate transformation, x = Aδ e, from orbital element differences to Cartesian coordinates are used [40]. A calibrated initial condition is calculated as follows.
24
69
72
1 T 1 T 1 T 1 v q˜ , q = q˜ 1q˜ + q˜ 1q˜ + q˜ K H B q˜ + q˜ T K H B q˜ 2 2 2 2 T 1 T 1 = −G q˜ − K H B q˜ 1q˜ + q˜ 1 −G q˜ − K H B q˜ (55) 2 2
˜
68 70
>0
which is a radially positive definite function. Taking the unbounded derivative of v q˜ , q˜ gives
˜
67
(52)
Due to the nature of eigenvalues then the system is stable. Therefore, the system is globally exponentially stable because the poles of the system will always have negative real parts. The total energy of the system given as
1 T 1 T v˜ q˜ , q˜ = q˜ 1 q˜ + q˜ K H B q˜ 2 2
(53)
Case 1: The equations of motion in radial, along-track and cross-track directions are integrated using ode45. The orbital elements of the chief and deputy are shown in Table 1 and the orbit propagation was performed for six orbits. The simulations are carried out using two different initial conditions; true and CW calibrated initial conditions.
121 122 123 124 125 126 127 128
(a) Uncalibrated Generating Solution In this simulation, the CW, HB, cubic and nonlinear models are propagated using true initial conditions, and the vector of constants c is also evaluated using true initial conditions. The simulation results are illustrated in
129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.7 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
7
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52 54 55 56
Table 1 Chief and Deputy orbital elements for a spacecraft in LEO orbit.
57
Orbital elements
Chief
Deputy
58
Semimajor axis, a (km) Eccentricity, e Inclination, i (deg) RAAN, (deg) Argument of Perigee, ω (deg) True Anomaly, f (deg)
7500 0 4 5 10 25
7500 0.0012 4 5.01 15 20
59 60 61 62 63 64 65 66
118
Fig. 3. Propagated solutions for Example 1 using uncalibrated generating solution.
53
Fig. 3. These results indicate that the HB model approximates the nonlinear dynamics more accurately than the CW model, at least
119
over some initial time interval. To emphasize this another plot of approximation errors, focusing on the first 1000 seconds is shown. To quantify the approximation errors, the following error metric is defined for the average error.
ε(t ) = ρˆ (t ) − ρ (t ), ε¯ =
1 t sim
t sim
ε(t )dt
120 121 122 123 124 125
(58)
126 127
0
128
where ρ (t ) is the nonlinear solution for the relative position, ρˆ (t ) is one of the approximate linear solutions, and t sim is the length of time that has been simulated. The average error for each linear solution is shown in Table 2.
129 130 131 132
JID:AESCTE
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.8 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
8
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114 115
49
Fig. 4. Propagated solutions for Example 1 using calibrated generating solution.
50
116 117
51 52 53 54
Table 2 Errors over 1000 seconds for Example 1 using uncalibrated generating solution.
Table 3 Errors over 1000 seconds for Example 1 using calibrated generating solution.
118 119 120
55
Model
Average error (km)
Model
Average error (km)
121
56
CW HB
0.325 0.040
CW HB
0.962 0.0484
122
57
125
59 60 61 62 63 64 65 66
123 124
58
(b) Calibrated Generating Solution In this simulation, the HB, cubic and nonlinear models are propagated using true initial conditions, while the vector of constants c is evaluated using calibrated initial conditions. For purposes of a fair or complete comparison, the CW solution is propagated using the calibrated initial condition. The simulation results are illustrated in Fig. 4. The values for ε¯ for each linear model
is shown in Table 3. These results illustrate the improvement in accuracy of the HB model when using the calibrated generating solution, compared with using the uncalibrated generating solution. When using the uncalibrated generating solution, the HB model exhibited erroneous drift. When using the calibrated generating solution, the HB model maintained better accuracy over a larger time span. However, averaged over the simulation interval, the HB
126 127 128 129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.9 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
9
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116 117
51 52
Fig. 5. Propagated solutions for Example 2 using uncalibrated generating solution.
54 55 56 57 58 59
Table 4 Errors over 1000 seconds for Example 2 using uncalibrated generating solution. Model CW HB
0.325 0.040
61 63 64 65 66
ing the calibrated generating solution is able to provide good long term error without having to introduce initial error, unlike the calibrated CW solution.
model using calibrated generating solution was less accurate than the calibrated CW solution. The tradeoff, though, is the initial error introduced in the calibrated CW solution. The HB model has no such initial error. This is illustrated in the initial approximation errors. The HB model us-
120 121 122 123
Average error (km)
60 62
118 119
53
Case 2: A second example is constructed with initial conditions defined identical to Table 1 except the deputy’s semi-major axis is changed to 7505 km. The CW, HB, cubic and nonlinear models are propagated using true initial conditions, and the vector of constants c is also evaluated using true initial conditions. Fig. 5 shows the propagated solutions using the uncalibrated generating solution. The HB method provides good accuracy initially, but the error appears to grow rapidly later in the simulation
124 125 126 127 128 129 130 131 132
JID:AESCTE
10
AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.10 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51 52 53 54 55 56 57
Table 5 Errors over 1000 seconds for Example 2 using calibrated generating solution.
59 61 62 63 64 65 66
(a) Calibrated generating solution The HB, cubic and nonlinear models are propagated using true initial conditions, while the vector of constants c is evaluated using calibrated initial conditions. The CW solution is propagated using the calibrated initial condition. Fig. 6 illustrates drifting orbit propagated solutions using calibrated generating solution and Table 5 shows errors over 1000 seconds.
118 119
interval. This may be related to how the HB model tries to approximate the nonlinear system as a linear, time varying system. To focus on the initial accuracy, the errors over the first 1000 seconds are shown in Table 4.
58 60
117
Fig. 6. Propagated solutions for Example 2 using calibrated generating solution.
Model
Average error (km)
CW HB
0.964 0.048
120 121 122 123 124 125 126
7. Conclusion
127 128
The harmonic-balance model of relative motion between chief spacecraft, in circular orbit, and deputy spacecraft has been presented. It results in an initial-condition dependent, linear model for the relative motion. The model is time-invariant. This paper has
129 130 131 132
JID:AESCTE AID:105761 /FLA
[m5G; v1.261; Prn:5/02/2020; 15:29] P.11 (1-11)
A.D. Ogundele / Aerospace Science and Technology ••• (••••) ••••••
1 2 3 4 5 6 7 8
shown that harmonic-balance model of spacecraft relative motion gives a more accurate linearized model because it captured the motion better than the conventional method, Clohessy-Wiltshire model. The results obtained from the numerical simulations show that harmonic-balance model has lesser error than the ClohessyWiltshire model of satellite relative motion. Also, from the stability analysis using total energy as a Lyapunov function, we obtained that the harmonic balance model is Lyapunov stable.
9 10
8. Acknowledgment
11 12 13 14
The author would like to acknowledge Dr. Andrew Sinclair of Air Force Research Lab and Dr. Subhash Sinha of Auburn University for explaining the concept.
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] G. Hill, Researches in lunar theory, Am. J. Math. 1 (1878) 5–26. [2] W. Clohessy, R. Wiltshire, Terminal guidance system for satellite rendezvous, J. Guid. Control Dyn. 27 (9) (1960) 653–658. [3] S. Lee, I. Hwang, Reachable set computation for spacecraft relative motion with energy-limited low-thrust, Aerosp. Sci. Technol. 77 (2018) 180–188. [4] V. Pescea, R. Opromollab, S. Sarnoc, M. Lavagnaa, M. Grassib, Autonomous relative navigation around uncooperative spacecraft based on a single camera, Aerosp. Sci. Technol. 84 (2019) 1070–1080. [5] R. Ambrose, I. Nesnas, F. Chandler, B. Allen, T. Fong, L. Matthies, R. Mueller, Nasa technology roadmaps: ta 4: Robotics and autonomous systems, Technical report, NASA. [6] Y. Wang, H. Ji, Integrated relative position and attitude control for spacecraft rendezvous with ISS and finite-time convergence, Aerosp. Sci. Technol. 85 (2019) 234–245. [7] J. Wertz, Orbit and Constellation Design and Management, Microcosm Press and Springer, New York, 2009. [8] M. Macdonald, V. Badescu, The International Handbook of Space Technology, Springer-Verlag, Heidelberg, New York, 2014. [9] H. Curtis, Orbital Mechanics for Engineering Students, 2nd edition, Elsevier Ltd, 2010. [10] H. Schaub, J. Junkins, Analytical Mechanics of Space Systems, 3rd edition, AIAA Education Series, 2014. [11] G. Gomez, M. Marcote, Lindstedt-Poincaré solutions of Clohessy-Wiltshire equations, in: AAS/AIAA Astrodynamics Specialists Conference, 2005, pp. 1–13, Paper AAS 05-359. [12] D. Lee, G. Vukovich, Kinematically coupled spacecraft relative motion without attitude synchronization assumption, Aerosp. Sci. Technol. 45 (2015) 316–323. [13] R. Battin, Astronautical Guidance, McGraw-Hill, New York, 1964. [14] R. Battin, An Introduction to the Mathematics and Methods of Astrodynamics, AIAA Education Series, New York, 1987.
11
[15] T. Carter, M. Humi, Fuel-optimal rendezvous near a point in general Keplerian orbit, J. Guid. Control Dyn. 10 (6) (1987) 567–573. [16] H. London, Second approximation to the solution of the rendezvous equations, AIAA J. 1 (1963) 1691–1693. [17] Y. Zhang, L. Yang, Y. Zhu, H. Huang, W. Cai, Modeling and analysis of dynamics for spacecraft relative motion actuated by inter-satellite non-contacting force, Aerosp. Sci. Technol. 43 (2015) 236–244. [18] J. Li, C. Gao, C. Li, W. Jing, A survey on moving mass control technology, Aerosp. Sci. Technol. 82–83 (2018) 594–606. [19] J. Sullivan, S. Grimberg, S. DAmico, Comprehensive survey and assessment of spacecraft relative motion dynamics models, J. Guid. Control Dyn. 40 (2017) 1837–1859. [20] I. Huntley, R. Johnson, Linear and Nonlinear Differential Equations, Ellis Horwood Series, 1983. [21] D. Jordan, P. Smith, Nonlinear Ordinary Differential Equations: An Introduction for Scientists and Engineers, Oxford University Press, 2011. [22] A. Nafeh, D. Mook, Nonlinear Oscillations, John Wiley and Sons, New York, 1979. [23] S. Yamgoue, On the harmonic balance with linearization for asymmetric single degree of freedom non-linear oscillators, J. Nonlinear Dyn. 69 (2012) 1051–1062. [24] A. Belendez, E. Gimeno, M.L. Alvarez, M. Yebra, D. Mendez, Analytical approximate solutions for conservative nonlinear oscillators by modified rational harmonic balance method, Int. J. Comput. Math. (2020). [25] Z. Guo, A. Leung, H. Yang, Iterative homotopy harmonic balancing approach for conservative oscillator with strong odd-nonlinearity, J. Appl. Math. Model. 35 (2011) 1717–1728. [26] M. Ghadimi, H. Kaliji, Application of the harmonic balance method on nonlinear equations, World Appl. Sci. 22 (4) (2013) 532–537. [27] E. Grebenikov, Encyclopedia of Mathematics, Kluwer Acad. Publrs, ISBN 1402006098, 2011. [28] A. Ogundele, Nonlinear dynamics and control of spacecraft relative motion, Ph.D. thesis, Auburn University, Auburn, Alabama, 2017. [29] A. Ogundele, A. Sinclair, S. Sinha, Developing harmonic balance model of spacecraft relative motion, in: AAS/AIAA Spaceflight Mechanics Meeting, Napa, CA. [30] K. Huseyin, Vibrations and Stability of Multiple Parameter Systems, Sijthoff and Noordhoff International Publishers, The Netherlands, 1978. [31] W. Thompson, G. Tait, A Treatise on Natural Philosophy, 1st edition, Cambridge University Press, Cambridge, England, 1879. [32] W. Haddad, V. Chellaboina, Nonlinear Dynamical System and Control: A Lyapunov Based Approach, Princeton University Press, 2008. [33] N. Butenin, Elements of the Theory of Nonlinear Oscillations, Blaidsdell Publishing Company, New York, 1965. [34] R. Shieh, Energy and variational principles for generalized (gyroscopic) conservative problems, Int. J. Non-Linear Mech. 55 (1971) 495–509. [35] K. Huseyin, R. Plaut, Transverse vibrations and stability of systems with gyroscopic forces, J. Struct. Mech. 3 (2) (1974) 163–177. [36] K. Huseyin, R. Plaut, Standard forms of the eigenvalue problems associated with gyroscopic systems, J. Sound Vib. 45 (1) (1976) 29–37. [37] D. Koditschek, The application of total energy as a Lyapunov function for mechanical control systems, Contemp. Math. 97 (1989) 131–157. [38] P. Hughes, Spacecraft Attitude Dynamics, Dover Publication, Inc., New York, 2004. [39] A. Sinclair, R. Sherrill, T. Lovell, Calibration of linearized solutions for spacecraft relative motion, J. Guid. Control Dyn. 37 (2014) 1362–1367. [40] K. Alfriend, H. Schaub, D. Gim, Gravitational perturbations, nonlinearity and circular orbit assumption effects on formation flying control strategies, in: AAS Rocky Mountain Guidance and Control Conference, vol. 104, 2000, pp. 139–158.
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65
131
66
132