Relative motion coupled control for formation flying spacecraft via convex optimization

Aerospace Science and Technology 14 (2010) 415–428

Contents lists available at ScienceDirect

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

Relative motion coupled control for formation flying spacecraft via convex optimization Yun-Hua Wu ∗,1 , Xi-Bin Cao 2 , Yan-Jun Xing 1 , Peng-Fei Zheng 1 , Shi-Jie Zhang 3 Research Center of Satellite Technology, Harbin Institute of Technology, 150080, Harbin, China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 10 October 2008 Accepted 13 April 2010 Available online 18 April 2010

Formation initialization control is of paramount importance for building a drift-free relative orbit, and is a challenging problem due to the coupled translational and rotational dynamics, e.g. the orientation of the thrust vector is constrained by the attitude and its angular velocity. We establish a nonlinear coupled dynamic model for formation flying spacecraft, and develop a relative orbit and attitude controller. The attitude angular velocity induced thrust vector constraint is then converted into thrust vector maneuverability constraint, which is nonconvex and cannot be implemented in the optimization framework directly. Thus the orbit and attitude controller can be designed separately. When designing the relative orbit controller, we use a relaxation method to convexify the nonconvex constraint and tailor the optimization problem to a semidefinite program, because of its low complexity and the existence of deterministic convergence properties. Then a variable structure attitude controller is used to track the optimized thrust direction. The validity of the proposed approach is demonstrated in a typical application of a dual-spacecraft formation initialization. © 2010 Elsevier Masson SAS. All rights reserved.

Keywords: Formation flying Coupled control Nonconvex constraints Convex optimization

1. Introduction Spacecraft Formation Flying (SFF) generally refers to a set of spatially distributed spacecraft flying in formation, capable of interacting and cooperating with each other to accomplish some space missions which are impossible for a monolithic spacecraft [18]. The U.S. Air Force and NASA are planning a number of formation missions to perform reconnaissance, virtual co-observing, stereo imaging, passive radiometry, optical interferometry, and terrain mapping [14,9]. These missions typically require an initialization control; which is defined as finding a fuel optimal trajectory that takes the deputy spacecraft with given initial states to prescribed final states with thrust magnitude and state constraints, to build a driftfree relative orbit. The spacecraft formation initialization control is challenging for small spacecraft equipped with thrusters that can only generate thrust along one axis in the body frame and reaction wheels that can only provide limited attitude maneuverability, and the six Degree-Of-Freedom (DOF) coupled translational and rotational dynamics due to subsystem inter-dependencies, actuator induced disturbances, and natural system coupling [13]. Subsystem inter-dependencies coupling include the dependence of thrust orientation for orbit control system on the spacecraft attitude and

* 1 2 3

Corresponding author. Tel.: +86 45186 402357x8507; fax: +86 45186 416457. E-mail address: [email protected] (Y.-H. Wu). Doctoral Candidate, Research Center of Satellite Technology. Professor, Research Center of Satellite Technology. Associate Professor, Research Center of Satellite Technology.

1270-9638/$ – see front matter doi:10.1016/j.ast.2010.04.005

© 2010

Elsevier Masson SAS. All rights reserved.

the thruster installation matrix. Furthermore, the spacecraft must satisfy attitude restrictions on the region of the sky where certain spacecraft instruments can point or restrictions on pointing towards other spacecraft for communication or plume avoidance [6]. Actuator induced disturbances include torques generated by thruster misalignment with the body axis. Natural coupling is apparent in such concepts as variation of magnetic and gravity gradient torques as a function of both attitude and orbit. Due to the nonconvexity of some constraints and the coupled dynamics, the coupled control problem is particularly difficult and is of paramount importance [20]. In recent years, much work has been concentrated on distributed spacecraft formation control. Due to the simplicity and the low requirements of control precision, a majority of these papers assume the spacecraft to be a point mass, hence the relative orbit and attitude control systems are designed separately, and the solutions for optimal trajectories are generally solved using 3 DOF relative orbital dynamics without any treatment of attitude induced constraints. The study of dynamics and control with coupled translational and rotational dynamics has received scant attention in the current literature except Refs. [14,20,19] and [16]. In Ref. [20], an output feedback tracking control law with a high-pass filter employed to estimate the relative velocities is designed based on mutually coupled translational and rotational dynamics of formation flying spacecraft. Some other methods, including sliding mode control [19] and state-dependent Riccati equation method [16], are used to control the spacecraft position and attitude to maneuver around tumbling space debris. Furthermore, an adaptive nonlinear

416

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

controller was developed with unknown mass and inertial parameters in Ref. [14]. However, all these coupled models are derived without considering the control constraints imposed by thrust direction and its rate of change, which is the major coupling for SFF, and these controllers are all designed with the assumption that the spacecraft is able to thrust in any direction. Neglecting the thruster vector constraints, such as actuator saturation and attitude maneuver capability, the trajectory might not be feasible to track. In this paper, our objective is to design a fuel optimal formation controller for the coupled system. The outline of this paper is as follows. In the first section, we present a detailed derivation of the nonlinear coupled dynamic equation with consideration of the thruster installation matrix. Next, the coupling imposed by the limited attitude angular velocity is then converted into a nonconvex thrust maneuverability constraint, thus the controller can be divided into independent relative orbital and attitude tracking controllers. To tailor the original problem to a convex optimization framework, especially the Semidefinite Program (SDP) framework, we present a general convexification procedure for handling the nonconvex constraints as a second order cone constraint according to a relaxation method. An available variable structure attitude tracking controller [3] is used to track the optimized thrust vector. Finally, we conclude our research with a numerical simulation example that demonstrates the validity of the developed framework. 2. Coupled dynamic equations Though the orbit and attitude dynamics for all spacecraft are coupled, the relative orbit and attitude control system are designed separately in traditional formation control problems, due to low requirements for control precision. For traditional orbit maneuver, the orbit control-firing period is very short relative to the orbit periods because of the large thrust provided by the propulsion system, and the effect of attitude control can be minimized or neglected. While for formation flying spacecraft equipped with low-level thrust propulsion systems, the orbit control can be regarded as a continuous low thrust control and thus the coupled influence cannot be ignored [4]. With the development and application of new space technology, low level thrust propulsion systems and Earth pointing observation missions, it necessitates the coupled modeling of translational and rotational dynamics, e.g., the dynamics of a solar sail propelled spacecraft is six DOF [10,8]. The relative orbit and attitude are also coupled for a tethered spacecraft system due to the tether constraints [15]. In this section, we will give a detailed derivation of the nonlinear coupled dynamics for SFF via the independent relative orbit and attitude dynamics. 2.1. Coordinate frame Earth Centered Inertial (ECI) frame O X I Y I Z I : This frame, shown in Fig. 1, is a nonaccelerating reference frame, with its origin at the center of the Earth, its X I axis pointing to the vernal, and its Z I axis pointing to the north and parallel to the rotation axis of the Earth. The Y I axis is given by the right-hand rule. Chief Orbit (CO) frame O c X o Y o Z o : The origin of this frame is located at the center of mass of the chief spacecraft, with its Z o axis pointing towards the Earth, and its Y o axis is parallel but opposite to the orbital angular momentum vector. The X o axis completes a right-hand coordinate frame with Y o and Z o axes. For a circular orbit, the X o axis will point in the direction of motion. Relative Motion (RM) frame: Its origin is at the chief spacecraft position and its orientation is given by the vector {x, y , z}. The x axis points anti-nadir, while the z axis is parallel to the orbit momentum vector in the chief orbit normal direction. The y axis then completes a right-handed set with the x and z axes, as shown in Fig. 1.

Fig. 1. ECI and RM frame.

Body Fixed (BF) frame O X b Y b Z b : This frame is centered on the center of mass of the spacecraft, with its axes fixed to the body and aligned along the principal axes of the spacecraft. 2.2. Nonlinear relative orbital dynamics The formation studied in this paper is a dual-spacecraft formation, in which one, denoted as the chief spacecraft, provides the basic reference trajectory, and the other, denoted as the deputy spacecraft, flies in the neighborhood of the chief spacecraft according to a desired relative trajectory. If we suppose that the chief spacecraft doesn’t generate any thrust during the formation maneuvering, the two spacecraft are identical and flying in a disturbance-free environment, the nonlinear relative orbital dynamics of the deputy spacecraft relative to the chief is governed by the following nonlinear equation [5] ◦◦

◦

r + 2C r + Q =

udRM

(1)

md

where the superscript RM denotes the vectors are expressed in the RM frame; md is the mass of the deputy spacecraft; r = [x, y , z]T is the relative position vector; the circle over a variable represents the time-derivative of the corresponding variable relative to the RM frame; u dRM is the orbit control vector for the deputy spacecraft, and can be given by



I RM RM RM udRM = M RM udI = udx , udy , udz

T

(2)

where the superscript I denotes the vectors are expressed in the I is the transformation matrix from the ECI frame ECI frame; M RM to the RM frame. The Coriolis-like matrix C and the nonlinear term Q are as

⎡

⎤

0 −1 0 C = ωo ⎣ 1 0 0 ⎦ 0 0 0



Q =

μ μ

R+x

 R + r 3 y

 R + r 3

−

(3) 1 R2

˙ o, − xωo2 − y ω

˙ o, μ − y ω + xω 2 o

z

 R + r 3

T (4)

where, ωo is the orbital angular velocity of the chief; R I is the position vector from the origin of the ECI frame to the center of mass of the chief; R is a scalar denotes the magnitude of R I ; μ = 3.986 × 105 km3 /s2 is the gravitational parameter of the Earth; the dots over a variable represent the time-derivative of the corresponding variable relative to the ECI frame.

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

417

2.3. Nonlinear attitude dynamics

3. Controller design

If the attitude is represented by quaternion q¯ = [qT , q4 ]T and angular velocity ω = [ω1 , ω2 , ω3 ]T , where q = [q1 , q2 , q3 ]T is the vector part and q4 is the scalar part, ω is the absolute angular velocity of the BF frame with respect to the ECI frame expressed in the BF frame. The quaternion satisfies the normalization constraint q¯ T q¯ = qT q + q24 = 1. The attitude kinematical equation of motion is given by

Resource optimization is vital for the lifespan of a formation mission, as in an outer space where refueling is not possible and it is necessary to reduce the overall weight of the spacecraft. As the result, low-thrust propulsion is becoming quite popular as it is characterized by high specific impulse, which allows reducing the amount of propellant needed for a mission [12]. The use of lowthrust continuous thrust propulsion leads to the necessity to find the optimal continuous control law. The formation optimal control problem can be stated as minimizing the objective function related to controls while satisfying a series of constraints. For Eq. (11), the optimization control problem can be described as

◦

q¯ =

1 2

Ω(ω)¯q

(5)

where ω is the angular velocity with respect to the RM frame, and satisfies ω = ω − M (¯q)ωoRM , M (¯q) is the transformation matrix from the RM frame to the BF frame; Ω(ω) is defined as



Ω(ω) =

−[ω×] ω −ωT 0



t f (6)

where the 3 × 3 matrix [ω×] is cross product operator and has the form

⎡

0

[ω×] = ⎣ ω3 −ω2

⎤ ω2 −ω1 ⎦

−ω3 0

ω1

(7)

0

The attitude dynamic equation of motion of a rigid spacecraft is then given by

˙ + [ω×] J ω = T + T d Jω

(8)

where, T are the control torques, T d are the disturbance torques, and J is the inertia matrix. 2.4. Coupled translational and rotational dynamics equation The coupled translational and rotational dynamics can be obtained from the nonlinear relative orbital and attitude dynamics derived above. We assume that each spacecraft is modeled as a rigid body with a thruster that provides body-fixed forces in one direction and torques about three mutually perpendicular axes in the BF frame. With the assumption that the thrust vector is aligned with the −x direction of the BF frame, Eq. (1) can be rewritten as

M T (¯q)udBF ◦ r + 2C r + Q  = md

◦◦

(9)

where udBF = [−1, 0, 0]T ud ; M T (¯q) is the transpose of M (¯q), and is given by

M T (¯q)

⎡

q24 + q21 − q22 − q23

= ⎣ −2(q4 q3 − q1 q2 )

2(q4 q3 + q1 q2 )

2(−q4 q2 + q1 q3 )

q24 − q21 + q22 − q23

2(q4 q1 + q2 q3 )

2(−q4 q1 + q2 q3 )

q24 − q21 − q22 + q23

2(q4 q2 + q1 q3 )

⎤T ⎦

◦

⎤

r ⎢r⎥ ⎢ ⎥

⎡

⎤

◦

⎡

r 03×1 ⎥ ⎢ C (¯q) ⎢ −2C r − Q  ⎥ ⎢ ⎢ 1 ⎣ q¯◦ ⎦ = ⎣ ⎦ + ⎣ 04×1 Ω(ω)¯q 2 03×1 − J −1 ω × J ω + J −1 T d ω˙ ◦◦

◦

◦

 f (r , q¯, ω) + Γ



ad T



(12)

0

where 0  α  1 is a scalar gain, t f is the terminal temporal node, subject to

0  ad  amax

   T (t )  T max

(13) (14)

and other state constraints, such as initial and terminal state constraints, etc. The control constraint (14) is equivalent to the attitude angular velocity constraint

  ω(t )  ωmax

(15)

Assuming that the attitude control is achieved by reaction wheel, and the orbit control uses low-level propulsion system, we can directly choose α = 0 for the optimization control problem of Eq. (12), for the propellant of orbit control is not refuelable on orbit. Eq. (11) is a 13-dimensional coupled dynamic equation, in ◦ which f (r , q¯, ω) is nonlinear, and Γ consists part of the system state; furthermore, the control sampling time of the relative orbit and attitude control subsystem are usually different. As a result, the optimal control problem of Eq. (12) is a difficult problem. In this paper, we convert the attitude angular velocity constraint of Eq. (15) into constraints of thrust direction and its rate of change. Then the resulting optimal control problem can be divided into relative orbit controller and attitude tracking controller and designed separately, as discussed in the following subsections. 3.1. SDP based relative orbit controller with thrust constraints

⎤

03×3   03×3 ⎥ ⎥ ad 04×3 ⎦ T J −1 (11)

where ad = ud /md is the thrust acceleration, and C (¯q) satisfies

C (¯q)

    α  T (t ) + (1 − α )ad dt

(10)

We choose the relative position and velocity, the quaternion, and the attitude angular velocity as the system state variable, and then the coupled dynamic equation can be written as

⎡

min

t f , T (·)

T    = − q24 + q21 − q22 − q23 −2(q4 q3 + q1 q2 ) −2(−q4 q2 + q1 q3 )

The spacecraft formation initialization control is challenging because of the constraints imposed by the thrust direction and its rate of change. Formation flying control problems with nonconvex constraints are often perceived to be computationally difficult and potentially in the class of NP-hard problems. In this section, we formulate the nonlinear and nonconvex thrust constraints and introduce a relaxation method to convexify the nonconvex constraints. The thruster for the relative orbit control is assumed to be fixed in the spacecraft BF frame, so the spacecraft has to alter its attitude to change the thrust direction. However, the achievable rotation rates are limited due to the maximum torques provided by the spacecraft attitude control system, which means that the thruster of the spacecraft cannot fire in any direction in a neighboring time interval. Formation flying spacecraft are generally small and usually have physical and operational constraints, an example of a physical constraint is the reduced space and mass allotted for the

418

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

orbital propulsion system. Many of these small spacecraft use low thrust propulsion systems, thus the available thrust magnitude is limited. The above constraints will be considered in designing the relative orbit controller, and are classified as maneuverability and magnitude constraints, respectively. The magnitude constraint is modeled as









U c = u (t ) ∈ R3 : 0  u (t )  u + , ∀t ∈ [0, t f ]

(16)

where R3 is the space of 3-dimensional vectors with real components, and u + ∈ R is the upper bound of thrust magnitude, with u + > 0. The initial and terminal state constraints are specified as ◦

◦

r (t 0 ) = r 0 ,

r (t 0 ) = r 0

r (t f ) ∈ r t f ,

r (t f ) ∈ r t f

◦

(17)

◦

(18) ◦

where t 0 is the initial temporal node; r t f and r t f are the relative position and velocity terminal constraint set, respectively. For the zero terminal constraint case, the terminal state constraint is ◦ an equality constraint r t f = r t f = 03×1 , while for the nonzero case, the terminal constraint can be described as a desired relative position terminal box

r t f ∈ C(δ)  {rt f ,i ∈ r t f : −δ  rt f ,i  δ, i = 1, 2, 3}

(19)

and a desired relative velocity terminal box ◦

◦

◦

(20)

Here we use terminal box rather than a generic norm constraint to represent the terminal constraints of Eqs. (19) and (20). The generic norm constraint, representing a constrained sphere, generally concerns the overall error, while the terminal box constraint deals with 3-axis error. Actually, both constraints can be used, but we are accustomed to choose the terminal box. Note that the error box of the terminal constraint should be chosen appropriately to save resources due to disturbances applied to the spacecraft or state estimation errors of the relative navigation system. If we assume that the maximum attitude angular velocity provided by the spacecraft attitude control system is ωmax with ωmax as its magnitude, and the body fixed thruster can only fire in one direction with respect to the BF frame. Thus, the change of the thrust vector during the time intervals t is limited by the following inequality

 cos θ(kt ) =

u T (kt )u ([k − 1]t )

u (kt ) · u ([k − 1]t )

k = 1, 2, . . .

min

t f ,u (·)

  u (t ) dt

0

subject to Eq. (7), and Eqs. (16)–(21).

min

Γ dt

t f ,u (·),Γ

(23)

0

where Γ satisfies

  u (t )  Γ

(24)

and

0  Γ  u+

(25)

For notational covariance, we introduce the following two variables

σ = Γ /md ,

a = udRM /md

(26)

For simplicity, we assume that the mass the spacecraft is invariant. Then, Eq. (6) can be rewritten as ◦◦

◦

r + 2C r + Q = a

(27)

Thus the original minimizing fuel problem is equivalent to minimize

t f (28)

0

with a(t )  σ (t ), and 0  σ (t )  amax , where amax denotes u + /md . The thrust maneuverability constraint (21) will be satisfied if the following constraint exists







a(kt )T a [k − 1]t  κσ (kt )σ [k − 1]t



(29)

where k = 1, 2, . . . , κ = cos(ωmax t ). By some algebraic operations on the inequality Eq. (29), we obtain



a(kt )T a(kt ) − a(kt )T a [k − 1]t



 T   + a [k − 1]t a [k − 1]t    σ (kt )σ (kt ) − κσ (kt )σ [k − 1]t     + σ [k − 1]t σ [k − 1]t

(30)

The above inequality can be written as



cos(ωmax t ) (21)

where u (kt ) is the control at the temporal node kt;  •  denotes the Euclidean norm on 3-dimensional Euclidean space, and θ(t ) is the net commanded angle of rotation for t ∈ [(k − 1)t , kt ]. The above inequality defines a nonlinear control constraint that cannot be directly implemented in an optimization control framework. Given the control constraints above, the minimum fuel formation optimization control problem considered in this paper can be described as

t f

t f

σ dt

◦

r t f ∈ C(ξ )  {rt f ,i ∈ r t f : −ξ  rt f ,i  ξ, i = 1, 2, 3}



With an additional slack variable Γ , the optimization problem (22) can be modified as [1]

a(kt ) a([k − 1]t )





I −0.5I



σ (kt ) σ ([k − 1]t )

−0.5I



I

aT (kt ) aT ([k − 1]t )

−0.5κ

1 −0.5κ



T T

σ (kt )

(31)

σ ([k − 1]t )

1

and can be simplified as

zk R zkT  sk P skT where



zk =

 (22)



R=

(32)



a(kt ) , a([k − 1]t ) I 3×3 −0.5I 3×3

 sk =

−0.5I 3×3 I 3×3



σ (kt ) σ ([k − 1]t ) 

,



P=

1 −0.5κ

−0.5κ



1

in which I 3×3 denotes the 3 × 3-dimensional identity matrix. Since |κ |  1, P = P T > 0, and R = R T > 0, therefore, we obtain

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

 1/2   1/2   R zk    P sk 

(33)

As λmin ( P 1/2 )sk    P 1/2 sk , consequently Eq. (33) will be satisfied if the following inequality is tenable

 1/2     R zk   λmin P 1/2 sk 

Case 1. sk  = 0, it is obvious that λmin ( P 1/2 )   P 1/2 sk . Case 2. sk  =  0, as P = P T > 0, there exists a matrix P 1/2 that 1/ 2 T 1/ 2 satisfies ( P ) P = P , and ( P 1/2 )T = P 1/2 for the matrix P in this paper. E.g.

P 1/2



=

1 4

√

√

4 − 2κ + 14 4 + 2κ √ √ −κ /( 4 − 2κ + 4 + 2κ )

√ √  −κ√ /( 4 − 2κ +√ 4 + 2κ ) 1 1 4 − 2κ + 4 4 + 2κ 4 (35)

Let λi ( P ) and e i be the eigenvalue and eigenvector of P . λi ( P ) and e i satisfy λ1 ( P ) > λ2 ( P ) > 0, P e i = λi ( P )e i and e Ti e j = 1 (i = j) or 0 (i = j) (i , j = 1, 2), respectively. Hence sk = a1 e 1 + a2 e 2 , for arbitrary sk ∈ R2 , with a21 + a22 = 0. Therefore, we have the square of the norm function  P 1/2 sk  as

 1/2 2  1/2 T 1/2    P sk  = P sk P sk = sT P 1/2 T P 1/2 sk k T   = P 1/2 P 1/2 sk , sk where ( P 1/2 )T P 1/2 sk = ( P 1/2 )T P 1/2

2

 P sk  sk  1/2

2

2

i =1 ai e i

=

  ( 2i =1 ai λi ( P )e i , 2i =1 ai e i ) =   ( 2i =1 ai e i , 2i =1 ai e i ) 2 2 =1 a i λ i ( P ) = i  λ2 ( P ) 2 2 i =1 a i

0

subject to Eqs. (24), (25), (27), and Eq. (41). SDP arises as convex relaxations of NP-hard quadratic optimization issues [2], and has been used to solve constrained attitude control problems [7]. SDP can be regarded as an extension of linear programming. Although SDP is much more general than linear programming, it is not much hard to solve. The advantages of SDP formulation include both the conceptual elegance and the efficient solvability via the available software [17]. In this section, we discretize the infinite-dimensional optimization problem (42) to finite-dimensional ones, which can be directly implemented in numerical algorithm. The discretization procedure is achieved by discretizing the time domain into equal time intervals and imposing the constraints at the temporal nodes. Sine the constraints are linear, quadratic, or second order cone constraints, the resulting problem is a finite-dimensional SDP that can be efficiently solved by SeDuMi [17] with a MATLAB interface [11]. The time interval [0, t f ] is divided into N equal subintervals [(i − 1)t , i t ], i = 1, . . . , N, where t is the length of the subinterval, thus the time nodes are given by

tk = kt ,

k = 0, 1 , . . . , N

(43)

The nonlinear differential Eq. (27) can be written in the following form

(36)

(44)

2

i =1 ai P e i

=

◦

◦ T

where X = [x, y , z, x, y , z] ; B = [03×3 , I 3×3 ] , a is the thrust acceleration, and Φ satisfies ◦

T

◦

Φ = [r , −2C r − Q ]T

(45)

The differential Eq. (44) can be easily converted into discrete form. The cost function (42) of the optimization problem can be approximated as follows

(37)

t f

σ dt ≈

min

t f ,a(·),σ (·)

      P 1/2 sk  λmin P 1/2 = λ2 P 1/2 = λ2 ( P )  sk 

(42)

◦

Using λi ( P ) to represent the eigenvalue and eigenvector of P 1/2 , we have (λi ( P 1/2 ))2 = λi ( P ), we have

(38)

min

t f ,a(·),σ (·)

0

N 

t σk

(46)

k =0

Thus, the original optimization problem is converted to discrete optimal control problem (46), subject to Eqs. (24), (25), (41), and Eq. (44). 3.2. Variable structure attitude tracking controller

As a result,

    λmin P 1/2 sk    P 1/2 sk 

(39)

Hence, the softening of inequality of Eq. (33) is proved. For optimal solution, we have σ ([k − 1]t ) ≈ σ (kt ), which leads to the following approximation

1   2

σ dt

min

t f ,a(·),σ (·)

◦

1/ 2

sk  ≈ √

t f

X = Φ + Ba

i =1 ai λi ( P )e i . We obtain



which is a second order cone constraint and can be implemented in the SDP framework directly. The minimum fuel formation initialization optimization problem can now be expressed as

(34)

where λmin ( P 1/2 ) is the minimum eigenvalue of P 1/2 . It is obvious that the inequality constraint of Eq. (33) will be satisfied if the inequality λmin ( P 1/2 )sk    P 1/2 sk  and Eq. (34) are tenable. Thus the inequality λmin ( P 1/2 )sk    P 1/2 sk  is the core between Eqs. (33) and (34), the detailed explanation is given as the following. We classify two cases to prove the above issue.

419





σ [k − 1]t + σ (kt )

 (40)

Then Eq. (31) can be approximately replaced by

 1/2  λmin ( P 1/2 )      R zk   σ [k − 1]t + σ (kt ) √ 2

The aim of the variable structure attitude tracking controller is to track the optimized thrust vector and its rate of change. The optimized thrust vector is described by quaternion qe , while its corresponding angular velocity ωe . Thus, the attitude tracking controller is given by [12]

T = [ω×] J ω + J

1 2



k sgn(δq4 ) Ξ T (¯q)Ξ (¯qe )ωe



˙ e − GΘ − Ξ (¯qe )Ξ (¯q)ω + ω T

(41)

with the sliding manifold

 (47)

420

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

s = (ω − ωe ) + k sgn(δq4 )δ q

(48)

where δq4 is the scalar part of the error quaternion δ q¯ defined by δ q¯ = q¯ ⊗ q¯ e−1 , in which the operator ⊗ denotes quaternion multiplication; k is a scalar gain; G is a 3 × 3-dimensional positive definite symmetric matrix; Ξ (¯qe ) is defined as Ξ (¯qe )  [qe,4 I 3×3 + [qe ×]T , −qe ]T , and Θ is a saturation function used to minimize chattering in the control torques and defined as

Θ k = sat(sk , εk ),

k = 1, 2, 3

(49)

in which Θ k and sk are the kth component of Θ and s respectively, and εk is a small positive quantity. The saturation function sat(sk , εk ) is defined by

sat(sk , εk ) =

⎧ ⎨1 ⎩

sk > εk

sk /εk

|sk |  εk , k = 1, 2, 3

−1

sk < −εk

(50)

The stability proof of the controller is omitted here, interested readers can refer to Ref. [12]. 4. Simulations, results, and discussions This section describes the application and verification of the method developed above for the SFF control problem. An example is presented to illustrate the initialization control for a dualspacecraft formation with coupled translational and rotational dynamics. The performance of the proposed method is investigated through numerical simulation. The physical constants, initial and terminal conditions used in our example are as follows. The inertia matrix along the principal axes is J = Diag([37.25, 18.32, 38.95]T ) kg · m2 . The initial attitude quaternion, angular velocity, relative position and velocity are given by

q¯ 0 = [0,

−0.8165,

−0.4082,

0.4082]T

ω0 = [0, 0, 0]T rad/s r 0 = [−655.2854, v 0 = [0.06,

1.07,

115.2300, 0.14]T m/s

−75.4500]T m

The desired relative translational trajectory is derived from the natural elliptical orbit of the chief and deputy around the Earth, and their orbit elements are chosen as in Table 1. The available maximum thruster acceleration amax , the maximum attitude angular velocity ωmax , the orbit control period T s (i.e. t), and the optimization steps N are chosen as amax = 0.002 m/s2 , which means the thruster can provide an F max = 0.4 N maximum thrust for a 200 kg satellite. If we suppose the specific impulse I sp of the thruster is 400 s, the orbit maneuver period is t = 1000 s, thus the estimated mass of fuel to be used is about m = t F max /( I sp g ) = 0.102 kg, which is really low (about 0.051%) compared to satellite mass, where g = 9.81 m/s2 is the gravitational acceleration. Hence, it can be assumed that the satellite mass is constant, ωmax = 0.8 deg/s, T s = 5 s, N = 200. Finally, the terminal box is chosen as δ = 5 cm, ξ = 0.4 mm/s for consideration of relative navigation errors. The attitude tracking control period T and other parameters are k = 0.06, G = 0.02I 3×3 , ε = 0.008, T = 0.5 s. The numerical simulation is executed as diagrammed in Fig. 2. The simulation results are shown from Fig. 3 to Fig. 10. Fig. 3 illustrates the 3-dimensional maneuver trajectory from the initial position marked as a star to the desired relative orbit described earlier, and the red line denotes the acceleration direction. The actual relative orbit represents the maneuver trajectory of the deputy spacecraft relative to the chief. The corresponding relative position and velocity tracking errors are shown from Figs. 4 and 5. Figs. 6 and 7 are the corresponding target quaternion and attitude angular velocity. Fig. 8 is the implemented thrust acceleration, with its magnitude and maneuverability constraints shown Table 1 Orbit elements for the chief and deputy spacecraft. Parameter and units

Chief

Deputy

Semimajor axis, km Eccentricity Inclination, deg Argument of perigee, deg RAAN, deg Mean anomaly, deg

6877.995832 0.001078 97.40617 90 85.7 0

6877.995815 0.00115 97.40617 90 85.7 0

Fig. 2. Coupled control scheme for formation flying spacecraft.

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

421

Fig. 3. 3-Dimensional sketch of the maneuver trajectory.

Fig. 4. Relative position tracking errors.

as Figs. 9 and 10, respectively. As Fig. 9 depicts, the thrust acceleration magnitude is strictly restricted under the magnitude bound. Fig. 10 shows the include angle of two thrust vectors in

a neighboring time interval, and it can be seen that the constraints imposed by the thrust direction maneuverability are satisfied.

422

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

Fig. 5. Relative velocity tracking errors.

Fig. 6. Target quaternion for tracking the thrust vector.

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

Fig. 7. Target attitude angular velocity for tracking the thrust vector.

Fig. 8. The implemented thrust acceleration.

423

424

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

Fig. 9. Magnitude of the implemented thrust acceleration.

Fig. 10. Include angle between thrust vectors in a neighboring time interval.

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

425

Fig. 11. Fuel consumption for different thrust acceleration (T s = 5 s, N = 200).

Fig. 12. Fuel consumption for different T s and N.

In order to analyze the relationship between fuel consumption and thrust acceleration bound, optimization steps and sampling time, we choose a different set of amax , T s and N, and give the simulation results from Fig. 11 to Fig. 14. Fig. 11 illustrates the fuel consumption curve vs. amax , which implies that the fuel consumption decreases as amax increases, but soon levels off after a

certain amax . Thus, a proper thruster can be chosen according to mission requirements. Fig. 12 is the sketch of the fuel consumption with different T s and N, the plot shows that the fuel consumption decreases with N or T s increasing. But the fuel consumption is invariable for a constant N × T s as shown in Fig. 13. Fig. 13 also implies that the fuel consumption is almost fixed for a speci-

426

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

Fig. 13. Fuel consumption for different N × T s .

Fig. 14. Include angle between thrust vectors for different N (T s = 5 s).

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

fied maneuver time, but decreases as the maneuver time increases. Thus, in order to reduce the computational load, we can choose a relatively small optimization steps and a relatively large time step.

427

Fig. 14 shows the include angle for different optimization steps. The plot indicates the include angle increases as N decreases, which means that the spacecraft has to fire agilely. As a result, the thrust direction maneuverability constraints may be

Fig. 15. Initial relative velocity error with different shooting angle.

Fig. 16. Fuel consumption for initial relative velocity with different shooting angle.

428

Y.-H. Wu et al. / Aerospace Science and Technology 14 (2010) 415–428

violated, e.g., the include angle violates the constraints with a small percentage when N = 140, which implies that the maneuver cannot be achieved under these constraints. To guarantee the optimization trajectory matches the constraints, a small margin should be added to N or T s before running the optimization algorithm. In order to analyze the effects of initial conditions, especially the initial relative velocity, on the fuel consumption, we perform simulations. We change the direction of the initial relative velocity, but the magnitude of relative velocity error and the initial relative position are chosen as previous examples. The initial relative position error r e0 and velocity error v e0 are defined as r e0 = r 0 − r d0 and v e0 = v 0 − v d0 , in which r d0 and v d0 are the desired relative position and velocity of the deputy spacecraft relative to the chief, respectively. Shooting Direction (SD) frame O X s Y s Z s is defined to describe the direction of the initial relative velocity error. The origin of the SD frame is at the initial relative position, with its X s axis pointing from the desired initial relative position to the actual initial position, and Z s axis is parallel to the direction of r e0 × v e0 . The Y s axis is given by the right-hand rule. We use shooting angle to describe the angle from Y s axis to the direction of the initial relative velocity error, seen in Fig. 15. The corresponding fuel consumption, shown as Fig. 16, is a sinusoid-like curve. It is obvious that the fuel consumption reaches maximum and minimum near 90◦ and 270◦ , respectively, which implies that the fuel consumption will reduce as the initial relative velocity direction points to the desired position. 5. Conclusions In this paper, we confine our attention to fuel optimal formation initialization control for a dual-spacecraft formation with coupled translational and rotational dynamics. The primary contributions of this work include: (1) A coupled relative dynamic equation was established; (2) The optimization problem for the coupled system was converted into two independent controllers, relative orbit controller and attitude tracking controller, by introducing a nonconvex thrust vector constraint; (3) The nonconvex constraint was convexified through a relaxation method, and the infinite-dimensional optimization problem was tailored to a finite-dimensional SDP optimization framework. A numerical simulation was presented to illustrate the efficacy of the proposed controller. Acknowledgements Research for this paper is supported by the National Natural Science Foundation of China, grant 60704020. The authors also would like to thank Dr. Dan Xue, Louis Scott Breger for several fruitful discussions on formation control, and Craig Mandsager for the spelling and grammar suggestions.

References [1] A.B. Acikmese, S.R. Ploen, A powered descent guidance algorithm for Mars pinpoint landing, in: AIAA Guidance, Navigation, and Control Conference and Exhibit, San Francisco, 2005, pp. 4515–4531. [2] S. Boyd, L. Vandenberghe, Semidefinite programming relaxations of non-convex problems in control and combinatorial optimizations, in: Communications, Computation, Control and Signal Processing: A Tribute to Tomas Kailath, Kluwer, 1997, pp. 279–288. [3] J.L. Crassidis, S.R. Vadali, F.L. Markley, Optimal variable-structure control tracking of spacecraft maneuvers, Journal of Guidance, Control and Dynamics 23 (3) (2000) 564–572. [4] A. Damiano, M. Bisten, Modeling attitude thrusters activity from attitude-orbit coupling effects, in: The 5th International Symposium on Spacecraft Operations, 1998. [5] M.S. de Queiroz, V. Kapila, Q. Yan, Adaptive nonlinear control of multiple spacecraft formation flying, Journal of Guidance, Control and Dynamics 23 (3) (2000) 385–390. [6] I. Garcia, J.P. How, Trajectory optimization for satellite reconfiguration maneuvers with position and attitude constraints, in: Proceedings of the American Control Conference, vol. 2, 2005, pp. 889–894. [7] Y. Kim, M. Mesbahi, G. Singh, F.Y. Hadaegh, On the constrained attitude control problem, in: AIAA Guidance, Navigation, and Control Conference, vol. 3, Providence, RI, United States, 2004, pp. 1862–1888. [8] D.A. Lawrence, S.W. Piggott, Integrated trajectory and attitude control for a four-vane solar sail, in: AIAA Guidance, Navigation, and Control Conference, 2005, pp. 2416–2431. [9] J. Leitner, F. Bauer, D. Folta, M. Moreau, R. Carpenter, J.P. How, Distributed spacecraft systems develop new GPS capabilities, in: GPS World: Formation Flying in Space, 2002, pp. 22–31. [10] M.E. Lisano, A practical six-degree-of-freedom solar sail dynamics model for optimizing solar sail trajectories with torque constraints, in: AIAA Guidance, Navigation, and Control Conference, 2004, pp. 828–838. [11] J. Löfberg, YALMIP: A toolbox for modeling and optimization in MATLAB, in: Proceedings of the IEEE International Symposium on Computer-Aided Control System Design, Taiwan, 2004, pp. 284–289. [12] Mauro Massari, Franco Bernelli-Zazzera, Optimization of low-thrust trajectories for formation flying with parallel multiple shooting, in: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Keystone, CO, 21–24 August 2006, pp. 1882–1897. [13] B.J. Naasz, M.M. Berry, H.-Y. Kim, C.D. Hall, Integrated orbit and attitude control for a nanosatellite with power constraints, Advances in the Astronautical Sciences 114 (2003) 1–18. [14] H. Pan, V. Kapila, Adaptive nonlinear control for spacecraft formation flying with coupled translational and attitude dynamics, in: Proceedings of the 40th IEEE Conference on Decision and Control, Orlando, FL, USA, 2001, pp. 2057– 2062. [15] S. Sgubini, P. Teofilatto, Attitude-orbit dynamical coupling in tethered systems, in: 6th International Conference on Dynamics and Control of Systems and Structures in Space, 2004. [16] D.T. Stansbery, J.R. Cloutier, Position and attitude control of a spacecraft using the state-dependent Riccati equation technique, in: Proceedings of the American Control Conference, Chicago, USA, 2000, pp. 1867–1871. [17] J.F. Sturm, Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones, Optimization Methods and Software 17 (6) (2002) 1105–1154. [18] C. Sultan, S. Seereeram, R.K. Mehra, F.Y. Hadaegh, Energy optimal reconfiguration for large scale formation flying, in: Proceedings of the 2004 American Control Conference, vol. 4, Boston, MA, 2004, pp. 2986–2991. [19] F. Terui, Position and attitude control of a spacecraft by sliding mode control, in: Proceedings of the American Control Conference, Philadelphia, PA, 1998, pp. 217–221. [20] H. Wong, H. Pan, V. Kapila, Output feedback control for spacecraft formation flying with coupled translation and attitude dynamics, in: Proceedings of the 2005 American Control Conference, 2005, pp. 2419–2426.