Relative orbit design and control of formation around displaced solar orbits

Aerospace Science and Technology 12 (2008) 195–201 www.elsevier.com/locate/aescte

Relative orbit design and control of formation around displaced solar orbits Shengping Gong ∗ , Hexi Baoyin, Junfeng Li Tsinghua University, Beijing, 100084, China Received 14 June 2006; received in revised form 25 May 2007; accepted 25 May 2007 Available online 2 June 2007

Abstract This paper investigates the relative orbit design and control methods for formation flying around displaced solar orbits. Firstly, the relative equations of motion are developed and linearized. Then, two types of orbit design and control methods are discussed based on the linearized equations. One method uses only position feedback to control the formation, in which the relative orbit can be controlled to be an ellipse. Further, the period of relative motion and direction of relative motion plane can be designed as mission requirements. The advantage of the method is that it employs an approximate system instead of a true one to design the parameters, which overcomes the trouble of solving a transcendental equation. The other method uses the full state feedback to control the formation, in which the relative orbit can be designed arbitrarily. Two control strategies, Linear Quadratic Regulation (LQR) and Input Feedback Linearization (IFL), are discussed for the latter. Finally, simulations are given to validate the two methods. © 2007 Elsevier Masson SAS. All rights reserved. Keywords: Solar sail; Displaced solar orbit; Formation flying; Orbit design

1. Introduction Solar sail has been investigated for several decades. Recently, several missions using solar sail have been proposed, such as new artificial L1 Lagrange point mission, which is used for NOAA/NASA to provide early warning of solar plasma storms, before they reach Earth [1,2]. Displaced solar orbit are another set of applications for solar sails. There are several prior references with regard to such orbits in early literatures. Recently, Forward considered a displaced solar sail north or south of the geostationary ring [3,4]. McInnes has done much work on this problem, in which large families of displaced orbits were discussed by considering the dynamics of a solar sail in a rotating frame, and the dynamics, stability and control of different families of displaced orbits of solar sails were investigated in detail [5–7]. Based on McInnes’s work, Molostov considered a more realistic solar sail model with non-perfect reflectivity and discussed the effect of finite absorption of the solar sail on the

* Corresponding author.

E-mail addresses: [email protected] (S. Gong), [email protected] (H. Baoyin), [email protected] (J. Li). 1270-9638/$ – see front matter © 2007 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.ast.2007.05.004

displaced orbits [8,9]. But the studies of the relative motion of solar sail around displaced orbit are rare in the literatures. This paper investigates the relative motion around displaced solar orbits. The paper is organized as follows: the second part introduces the displaced solar orbit briefly. The third part derives the relative equations of motion around a displaced orbit and the equations are linearized in the vicinity of the displaced orbit. The last two parts discuss the orbit design and control methods. The fourth part investigates the orbit design method when only the position feedback is employed to control the formation. An approximate system is introduced to design the relative orbit, which overcomes the trouble of solving a transcendental equation. The last part discusses the full state feedback control methods. LQR and IFL are used to track the reference relative orbit, and simulations show that they yield similar results. 2. Displaced solar orbit The displaced solar orbit is a circular Sun-centered orbit displaced above the ecliptic plane by directing a component of the solar radiation pressure force normal to the ecliptic. The sail orientation is defined by the vector n, fixed in the rotating

196

S. Gong et al. / Aerospace Science and Technology 12 (2008) 195–201

frame and sail performance is characterized by the sail lightness number β. To keep equilibrium in the rotating frame, the sail pitch angle α and sail lightness number β for a displaced Sun-centered orbit of radius ρ, displacement z and orbital angular velocity ω can be given by [6] tan α =

z ˜ 2 ρ (ω/ω) (z/ρ)2 + 1 − (ω/ω) ˜ 2

  2 1/2 {(z/ρ)2 + [1 − (ω/ω) ˜ 2 ]2 }3/2 z β = 1+ ρ [(z/ρ)2 + 1 − (ω/ω) ˜ 2 ]2

(5)

The transformation of the acceleration between the inertial and rotational frames can be given as

(1)

d 2δ ¨ = δ + 2ω × δ˙ + ω × (ω × δ) + ω˙ × δ dt 2

(2)

where d/dt denotes derivative with respect to the rotating frame. In this case, ω˙ = 0. Thus, the equation of motion in the rotating frame can be written as

where ω˜ 2 = μ/r13 .

(6)

δ¨ + 2ω × δ˙ + ω × (ω × δ) = f (r 1 ) − f (r 2 ) + g(r 2 , n2 )

3. Relative motion

− g(r 1 , n1 ).

In this paper, a two-sail formation is considered in which the leader revolves on a displaced solar orbit and the follower is on its nearby orbit. In this formulation, the sail’s state vector is defined in terms of a set of Cartesian coordinates relative to the rotating frame. The origin of this frame is the mass center of the Sun. In the rotating frame, the x-axis is in the ecliptic plane and rotates around the Sun with an angular velocity ω, the zaxis is in the direction of the angular velocity; and the y-axis forms a right-handed triad. Another rotating frame, labeled as the leader frame oξ ςη, has all its axes parallel to the axes of the rotating frame Oxyz and has its origin at the mass center of the leader. In the frame oξ ςη, the position vector of the follower can be written as δ = [ x2 − x1 y2 − y1 z2 − z1 ]T = [ ξ ς η ]T , as shown in Fig. 1. The equations of motion of the sails in an inertial frame can be written as d 2r 1 = −f (r 1 ) + g(r 1 , n1 ), dt 2 d 2r 2 = −f (r 2 ) + g(r 2 , n2 ) dt 2

d 2δ d 2r 2 d 2r 1 = − dt 2 dt 2 dt 2 = f (r 1 ) − f (r 2 ) + g(r 2 , n2 ) − g(r 1 , n1 ).

(3) (4)

where f (r) = (μ/r 3 )r is the solar gravitational force, g(r, n) = β(n · r)2 (μ/r 4 )n is the solar radiation pressure force, and μ is the solar gravitational constant. Then, the relative equation of motion in the inertial frame can be obtained by subtracting two equations above, yielding

(7)

The relative distance between sails |r 2 − r 1 | is usually very small compared with r1 or r2 . Therefore, f (r 2 ) can be expanded around r 1 . To maintain the follower in the vicility of the leader, the sail lightness number β2 should be very close to β1 and so should the normal vector n2 . Therefore, g(r 2 , n2 , β2 ) can also be expanded around (r 1 , n1 , β1 ). Then, the first order approximation of the relative equation of motion can be determined as δ¨ + 2J δ˙ + Mδ = N δn + P δβ

(8)

where r 1 = [ ρ 0 z ]T , n1 = [ cos γ1 0 sin γ1 ]T , in which γ is the sail elevation angle [6] measured from the ecliptic plane in the plane spanned by r and ω. The coefficient matrices are given by ⎡ ⎤ 0 −ω 0 J = ⎣ ω 0 0 ⎦, 0 0 0  2β1 μ(n1 · r 1 ) 4β1 μ(n1 · r 1 )2 M =− n1 nT1 − n1 r T1 4 r1 r16  μ 3μ − 3 I3×3 + 5 r 1 r T1 , r1 r1 N=

β1 μ(n1 · r 1 )2 2β1 μ(n1 · r 1 ) I3×3 + n1 r T1 , 4 r1 r14

P=

μ(n1 · r 1 )2 n1 . r14

δn can be further linearized by introducing two attitude angles, θ and φ, where θ is the angle between the normal and the xz plane, and φ is the angle between the x-axis and the projection of the normal in the xz plane. Then, the normal vector of the solar sail can be given by n = [ cos θ cos φ

Fig. 1. Coordinates and angles definition.

sin θ

cos θ sin φ ]T .

(9)

The control variables are defined by u = [ δφ δθ δβ ]T , δφ = φ2 − φ1 , δθ = θ2 − θ1 = θ2 and δβ = β2 − β1 . To maintain the follower in the vicinity of the displaced solar orbit, δφ and δθ should be small. Therefore, n2 can be expanded around (φ1 , θ1 ) and the first order approximation can be obtained as

S. Gong et al. / Aerospace Science and Technology 12 (2008) 195–201

  ∂n

δφ δn = n2 − n1 =

∂(φ, θ ) (φ1 ,θ1 ) δθ ⎡ ⎤ − sin φ1 0   δφ 0 1⎦ =⎣ . δθ cos φ1 0

be zero by selecting proper initial values. Then, Eq. (16) can be rewritten as (10)

The controlled relative equation of motion is obtained by substituting Eq. (10) into Eq. (8), yielding δ¨ + 2J δ˙ + Mδ = Bu.

(11)

The coefficient matrix is given by

  ∂n

P . B= N ∂(φ, θ ) (φ1 ,θ1 ) 4. Orbit design and control using position feedback Since the relative velocity is difficult to measure in engineering, only position feedback is utilized to control the formation in this section. The period of relative motion and the direction of relative motion plane are two important elements for formation flying. An orbit design method is used to design the two elements as the mission requirements. A generic position feedback controlled system can be given by δ¨ + 2J δ˙ − V δ = 0.

(12)

Because the damping coefficient matrix is antisymmetric, it can be proved that if the matrix V is symmetric negative definite, the whole linear system is Lyapunov stable, but not asymptotically stable [10,11]. For formation flying, attention is focused on the stable solutions, because the asymptotically stable solutions mean that the follower will impact with the leader. Therefore, the stable system can be obtained by designing V to be symmetric negative definite through position feedback. Then, the position feedback matrix Vr can be obtained as Vr = M + V .

(13)

The control law can be given by u = B −1 Vr δ.

(14)

The characteristic equation of the controlled system is given by  2 (15) λ I + 2J λ − V υ = 0 where I is a unit matrix, and λ is an eigenvalue, υ = υ r + υ i j is √ the corresponding eigenvector and j = −1. The system is stable, but not asymptotically, if V is symmetric negative definite [10]. Therefore, all the eigenvalues are pure imaginary or zero, and they can be written in the form of ±ωi j (i = 1, . . . , 3). Then, the solution of the controlled system can be given by δ=

3

      Ai cos ωi t + ψi0 υ ri − Ai sin ωi t + ψi0 υ ii

197

(16)

i=1

where Ai is the amplitude, and ψi0 is the initial phase, which is assumed to be zero hereafter. If all the three modes are included in the relative motion, the orbit is not periodic because ωi are not commensurable numbers. To achieve a periodic relative motion with angular velocity ωi . Two of the three amplitudes can

δ = Ai cos(ωi t)υ ri − Ai sin(ωi t)υ ii .

(17)

The relative motion described by Eq. (17) is restricted in the plane spanned by υ ri and υ ii , and the relative orbit is elliptical [10]. The parameters determining the shape of the ellipse can be given by     r2 (|υ ri |2 − |υ ii |2 )2  r i 2  |υ i | + |υ ii |2 + + υi · υi , a = Ai 2 4 (18)     r2 (|υ ri |2 − |υ ii |2 )2  r i 2  |υ i | + |υ ii |2 − + υi · υi , b = Ai 2 4 (19) √ a 2 − b2 . (20) e= a The normal of the relative motion plane is determined by the eigenvector, given by υ i × υ ri δ × δ˙ hˆ = . = ii ˙ |υ i × υ ri | |δ × δ|

(21)

The period of the relative motion is determined by the eigenvalue. The procedure of orbit design can be described as: first, the eigenvalue and eigenvector are obtained according to the period and the direction of relative motion plane required by the mission. Then, the position feedback matrix Vr is designed to guarantee that the controlled system has the corresponding eigenvalue and eigenvector. An eigenvalue of the controlled system is assumed to be ωc j and the corresponding eigenvector is υ r + υ i j . Let υ r = [ a1 a2 a3 ]T , and υ i = [ b1 b2 b3 ]T . Substitute them into Eq. (15) to obtain the relations given by ⎡ ⎤⎡ ⎤ 2ωc ωj + Vxy Vxz ωc2 + Vxx a1 + b1 j ⎢ ⎥ ωc2 + Vyy Vyz ⎦ ⎣ a2 + b2 j ⎦ ⎣ −2ωc2 ωj + Vxy a3 + b3 j Vxz Vyz ωc2 + Vzz = 0. (22) It can be rearranged as ⎡ ⎤⎡ ⎤ Vxx a1 a 2 a 3 0 0 0 ⎢ b1 b2 b3 0 0 0 ⎥ ⎢ Vxy ⎥ ⎢ ⎥⎢ ⎥ ⎢ 0 a1 0 a2 a3 0 ⎥ ⎢ Vxz ⎥ ⎢ ⎥⎢ ⎥ ⎢ 0 b1 0 b2 b3 0 ⎥ ⎢ Vyy ⎥ ⎢ ⎥⎢ ⎥ ⎣ 0 0 a1 0 a2 a3 ⎦ ⎣ Vyz ⎦ 0 0 b1 0 b2 b3 Vzz ⎡ ⎤ ⎡ ⎤ a1 −b2 ⎢ b1 ⎥ ⎢ a2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ b1 ⎥ 2 ⎢ a2 ⎥ ⎢ ⎥. = −ωc ⎢ ⎥ − 2ωc ω ⎢ ⎥ ⎢ b2 ⎥ ⎢ −a1 ⎥ ⎣ a3 ⎦ ⎣ 0 ⎦ b3 0 The last two rows of Eq. (23) are given by

(23)

198

  a1 Vxz + a2 Vyz + a3 ωc2 + Vzz = 0,   b1 Vxz + b2 Vyz + b3 ωc2 + Vzz = 0.

S. Gong et al. / Aerospace Science and Technology 12 (2008) 195–201

(24) (25)

From the equations above, it can be concluded that [Vxz Vyz ωc2 + Vzz ]T is perpendicular to both υ ri and υ ii , therefore, the normal of the relative motion plane can be obtained as  T  T

hˆ = Vxz Vyz ωi2 + Vzz / Vxz Vyz ωi2 + Vzz . (26) Eq. (23) is a transcendental equation of V because the right hand terms include the eigenvalue of V . Therefore, it is difficult to obtain the solution, even numerically. To avoid solving the transcendental equation, a simplified approximate system is introduced to design the position feedback matrix. The approximate system that ignores the Coriolis force is given by δ¨ − V δ = 0.

(27)

The angular velocity of rotating frame is smaller compared with the relative angular velocity, and the Coriolis force induced by the rotation of the frame will change the eigenvalues of the approximate system slightly. Therefore, the design results based on the approximate system are slightly different from the results of the true system, so the design method based on the approximate system could be improved because of the difference. In the subsequent analysis, two angles are employed to describe the direction of the relative orbit plane. Inclination i is the angle between the normal of the relative plane and z axis. Clock angle Ω is the angle between the x axis and the projection of the normal in the xy plane. Assume the design parameters are based on id , Ωd , and the relative angular velocity ωd , where the subscript d indicates the design parameters. 4.1. Design potential matrix based on the approximate system When the control variables u = 0, the potential matrix can be written as ⎡ ⎤ −M11 0 −M13 −M22 0 ⎦. V = −M = ⎣ 0 (28) 0 −M33 −M31 It can be concluded from Eqs. (26) and (28) that the normal always stays in the xz plane with u = 0. To guarantee that the normal can be arbitrary, the relations, Vxz = 0 and Vyz = 0, should be satisfied. Namely, the last row of the controlled potential matrix should be nonzero. To keep the sparseness of the matrix, other zero elements are kept. Here we design the potential matrix V to be symmetric, and then the controlled matrix can be written as ⎤ ⎡ Vxx 0 Vxz Vyy Vyz ⎦ . (29) V =⎣ 0 Vxz Vyz Vzz The eigenvalue of the approximate system is designed to be ωd j , therefore, the characteristic equation of Eq. (27) is written as

2

ω I + V = 0. (30) d

Eq. (30) can be expanded as     2 ωd + Vxx ωd2 + Vyy ωd2 + Vzz  2   2  2 2 ωd + Vyy − Vyz ωd + Vxx = 0. − Vxz

(31)

It can be found that one special solution to Eq. (31) can be given by Vxx = −ωd2 ,

(32)

Vyy = −ωd2 .

(33)

Eqs. (24) and (25) are still valid for the approximate system. According to the definition of inclination i, clock angle Ω and the normal of the relative motion plane hˆ given by Eq. (26), V can be designed as ⎧ ⎨ Vxz = k sin id cos Ωd , Vyz = k sin id sin Ωd , (34) ⎩ V = k cos i − ω2 zz d d where k is a constant. To guarantee the stability of controlled system, V is designed to be negative definite, and the necessary condition is given by |V | < 0.

(35)

Eq. (34) is substituted into Eq. (35) to obtain 2 2 − Vyz > 0. Vxx Vzz − Vxz

The solution is found to be  1 + 3 sin2 id − cos id 0 < k < |Vxx | . 2 sin2 id

(36)

(37)

Now, all the elements of the controlled potential matrix V have been determined. Therefore, the position feedback matrix is obtained. However, the eigenvalue and inclination of the true system are different from the design parameters. Therefore, the relation between the true system and the approximate system should be built up. 4.2. Map the true system and the approximate system The eigenvalue of the true system is ωt j instead of ωd j due to the Coriolis force, the inclination and clock angle of the true system are given by  tan Ωt = Vyz /Vxz , (38) cos it = (Vzz + ωt2 )/k. It can be found that Ωt = Ωd but it = id because Ωt is independent of the eigenvalue of the system and uniquely determined by the potential matrix. The inclination it , however, depends on both the potential matrix and eigenvalue of the system. Therefore, the angular velocity and inclination of the true system are different from the design parameters, and a numerical method is used to establish the relation between the true and approximate systems. The parameter pair (ω, i) is considered and a map between (ωd , id ) and (ωt , it ) is expressed as  (39) (ωd , id ) = F (ωt , it ) .

S. Gong et al. / Aerospace Science and Technology 12 (2008) 195–201

199

Fig. 2. (ωd , id ) pair vs. (ωt , it ) pair.

Fig. 3. Orbit design and control process.

Given the design pair (ωd , id ), the potential matrix V can be obtained based on the approximate system. Then, the Coriolis force is added to the system and the true pair (ωt , it ) can be calculated. Therefore, the map can be built up. Fig. 2 gives the results of (ωd , id ) pair vs. (ωt , it ) pair. Now, the design and control process can be described in Fig. 3. If a desired pair (ωt , it ) is given, as discussed above, the corresponding design pair (ωd , id ) can be obtained by interpolating the calculated data. Then, the controlled potential matrix V is designed based on the approximate system, and the position feedback matrix is obtained. Finally, the attitude and lightness number of the follower sail are calculated according to the control force.

Fig. 4. Relative orbit and control variables.

is 100 times of that of Earth’s; the desired inclination and clock angle are both π/4. Therefore, the desired pair (ωt , it ) is (1.99 × 10−6 , 0.7854). The corresponding design pair (ωd , id ) is found to be (2.06 × 10−6 , 0.7367), and the potential matrix V can be designed based on the method above. Fig. 4 gives the relative orbit and control law. The simulation shows that the designed results correspond with the desirable results very well. 5. Orbit design and control using position and velocity feedback

4.3. Numerical simulation In the following example, the integrations are based on the nonlinear equation of relative motion. The leader is on an Earth-synchronous displaced solar orbit with ρ = 0.98 AU and z0 = 4 × 105 km. The desired relative angular velocity

The difficulty in orbit design using position feedback becomes a trivial when the relative velocity is used to control the formation. All the typical control methods can be used to control the formation. In this part, two methods are discussed, Linear Quadratic Regulation (LQR) and Input Feedback Lin-

200

S. Gong et al. / Aerospace Science and Technology 12 (2008) 195–201

earization (IFL). They can be used to track the desired relative orbit. LQR controller is based on the linear model and the IFL is based on the nonlinear model. Therefore, results of the same relative orbit controlled by the two methods can be used to validate the linear model. The IFL controller will be not discussed here, and the simulation result is given. 5.1. LQR The first-order relative equation of motion is written as   0 ˙ Xr = AXr + u (40) B where Xr = [ δ T δ˙ T ]T . B and u are given in Eq. (11). It can be proved that the system is controllable. Consider the general quadratic cost function 1 J= 2

∞  (Xr − Xref )T Q(Xr − Xref ) + uT Ru dt

6. Conclusions (41)

0

where the matrices Q and R represent the weight of the state error and the control input, respectively, both of which are symmetric positive definite. Therefore, what is required is to find a linear state feedback by solving a Riccati matrix algebraic equation to get u = −K(Xr − Xref )

control inputs. Fig. 5A shows that the controlled relative orbits for two different methods are not completely consistent with each other. The maximum distance difference for the orbit is about several meters, which is small for a 10 km formation flying. Since the IFL controller considers the nonlinear system, the corresponding controlled orbit will be consistent well with the reference orbit. However, the LQR controller is based on the linear system, and the corresponding controlled orbit will be influenced by the nonlinear terms that are ignored by controller. Therefore, it can be concluded that the difference between the two controllers is introduced by the nonlinear terms ignored by the LQR controller. According to the time-history of the control variables, the formation requires a maximum control angle variation less than 1 degree so that the validity of the linearization can be guaranteed.

(42)

where Xref is the reference relative orbit. 5.2. Numerical simulation The reference relative orbit of the follower sail is a 10 km radius circular orbit in the ξ ς plane with its center at the origin of the frame oξ ςη. The period of the relative motion is one percent of Earth’s orbital period. The follower is assumed to be on the reference orbit at the initial time. Where Q = I6×6 and R = 1 × 106 I3×3 are set for LQR, the results obtained for the two control methods are shown in Fig. 5, where Fig. 5A illustrates the relative orbit and Figs. 5B, 5C and 5D show the

Based on the linearized equation of relative motion around a displaced orbit, the paper investigates the relative orbit design and control methods for the formation. A transcendental equation must be solved for the design method using only position feedback to control the formation. To avoid solving this equation, an approximate system is proposed to design the relative orbit. Then, a map between the approximate system and the true system is built up. Therefore, the design parameters referred to the true system can be transformed into the approximate system, on which all the design is based. The numerical examples show that the method can design the desirable relative orbit efficiently. The method utilizing full state feedback can track an arbitrary reference relative orbit. Therefore, the orbit design and orbit control are equivalent for this method. Two typical controllers, LQR and IFL, are employed to control the formation. The simulation results show that they yield similar results. Acknowledgements This work was supported by the National Natural Science Foundation of China (Grants No. 10602027 and No. 10672084). References

Fig. 5. Formation control with LQR and IFL.

[1] J.Y. Prado, A. Perret, G. Pignolet, I. Dandouras, Using a solar sail for a plasma storm early warning system, International Academy of Astronautics, Paper 96-IAA. 3.3.06, Oct. 1996. [2] J.L. West, NOAA/DOD/NASA Geostorm Warning Mission, Jet Propulsion Lab, California Institute of Technology, JPL D-13986, Pasadena, CA, Oct. 1996. [3] R.L. Forward, Light-levitated geostationary cylindrical orbits, Journal of the Astronautical Sciences 29 (29) (1981) 73–80. [4] R.L. Forward, Light-levitated geostationary cylindrical orbits using perforated light sails, Journal of the Astronautical Sciences 32 (4) (1984) 221–226. [5] C.R. McInnes, F.L. Simmons, Solar sail Halo orbits I: Heliocentric case, Journal of Spacecraft and Rockets 29 (4) (1992) 466–471. [6] C.R. McInnes, Passive control of displaced solar sail orbits, Journal of Guidance, Control, and Dynamics 21 (6) (1998) 975–982.

S. Gong et al. / Aerospace Science and Technology 12 (2008) 195–201

[7] C.R. McInnes, Solar sail mission applications for non-Keplerian orbits, Acta Astronautica 45 (1999) 567–575. [8] A.A. Molostov, A.A. Shvartsburg, Heliocentric Halos for a solar sail with absorption, Soviet Physics Doklady 37 (3) (1992) 149–152. [9] A.A. Molostov, A.A. Shvartsburg, Heliocentric synchronous halos for a solar sail with absorption, Soviet Physics Doklady 37 (4) (1992) 195–197.

201

[10] F.Y. Hsiao, D.J. Scheeres, The dynamics of formation flight about a stable trajectory, Journal of the Astronautical Sciences 50 (3) (2002) 269–287. [11] F.Y. Hsiao, D.J. Scheeres, Design of spacecraft formation orbits relative to a stabilized trajectory, Journal of Guidance, Control and Dynamics 28 (4) (2005) 782–794.