- Email: [email protected]
Switch programming of reflectivity control devices for the coupled dynamics of a solar sail
Available online at www.sciencedirect.com
ScienceDirect Advances in Space Research 57 (2016) 1147–1158 www.elsevier.com/locate/asr
Switch programming of reflectivity control devices for the coupled dynamics of a solar sail Tianjian Hu a,b,⇑, Shengping Gong a, Junshan Mu a,c, Junfeng Li a, Tianshu Wang a, Weiping Qian b b
a School of Aerospace Engineering, Tsinghua University, Beijing 100084, China Beijing Institute of Tracking and Telecommunication Technology, Beijing 100094, China c China Satellite Maritime Tracking and Control Department, Jiangyin 214431, China
Received 15 September 2015; received in revised form 16 December 2015; accepted 17 December 2015 Available online 25 December 2015
Abstract As demonstrated in the Interplanetary Kite-craft Accelerated by Radiation Of the Sun (IKAROS), reflectivity control devices (RCDs) are switched on or off independently with each other, which has nevertheless been ignored by many previous researches. This paper emphasizes the discrete property of RCDs, and aims to obtain an appropriate switch law of RCDs for a rigid spinning solar sail. First, the coupled attitude–orbit dynamics is derived from the basic solar force and torque model into an underdetermined linear system with a binary set constraint. Subsequently, the coupled dynamics is reformulated into a constrained quadratic programming and a basic gradient projection method is designed to search for the optimal solution. Finally, a circular sail flying in the Venus rendezvous mission demonstrates the model and method numerically, which illustrates approximately 103 km terminal position error and 0.5 m/s terminal velocity error as 80 independent RCDs are switched on or off appropriately. Ó 2015 COSPAR. Published by Elsevier Ltd. All rights reserved.
Keywords: Solar sail; Coupled attitude–orbit dynamics; Reflectivity control devices; Gradient projection method
1. Introduction Due to its capability of generating propulsion without fuel consumption, the solar sail has been envisioned to enable or enhance a wide range of mission applications throughout the solar system (Macdonald and McInnes, 2011). Some previous literature (Gong and Li, 2014; Zeng et al., 2014) mainly focused on the orbit design of a solar sail and assumed that the attitude can change
⇑ Corresponding author at: School of Aerospace Engineering, Tsinghua University, Beijing 100084, China. E-mail addresses: [email protected] (T. Hu), [email protected] (S. Gong), [email protected] (J. Mu), [email protected]. cn (J. Li), [email protected] (T. Wang), qianweipingbittt@sina. com (W. Qian).
http://dx.doi.org/10.1016/j.asr.2015.12.029 0273-1177/Ó 2015 COSPAR. Published by Elsevier Ltd. All rights reserved.
instantaneously. However, the propulsion exerted by solar radiation pressure (SRP) is related to the angle between the sunlight and the normal vector of the sail surface, so the orbit dynamics and attitude dynamics are not completely independent of each other (Gong et al. 2009; Zhang and Wang, 2013). Thus, neither the attitude maneuver process nor the attitude–orbit coupled effects can be neglected. Starting from the above idea, Gong et al. (2009) developed the coupled attitude–orbit dynamics for displaced solar orbits and discussed the stability of its motion. Zhang and Wang (2013) considered the flexibility of the sail membrane and derived a coupled dynamical equation for a flexible displaced solar orbiter. Nevertheless, neither of these authors discussed the physical realization of the desired trajectories in both orbit and attitude. An interesting implementation of reflectivity control devices (RCDs)
1148
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
on the attitude adjustment of a solar sail was demonstrated by the Japanese sail ‘‘interplanetary kitecraft accelerated by radiation of the sun” (IKAROS) (Tsuda et al., 2011a,b, 2013). RCD is a thin film attached to the surface of the sail that can switch its reflectivity between the specular and diffusive states versus electricity power being switched on or off. Different reflectivity states can result in different SRP. By adjusting the total amount and distribution of switched-on RCDs, the propulsion and torque acting on the sail can be controlled to follow the desired orbit and attitude. Compared to other methods for attitude control utilizing the gimbaled control boom (Wie 2004a,b), sail control vanes (Wie 2004a,b), or sliding masses (Scholz et al., 2011), the application of RCDs is obviously propellantless, much simpler in structure and has lower cost (Tsuda et al., 2011a,b). Such benefits have prompted increasing research interest in the utilization of RCDs in solar sail missions. Mu et al. (2013, 2015a,b) extended the coupled attitude–orbit dynamics to formation flying for a GeoSail mission. The reflectivity modulation ratio and its distribution on the sail are used as continuous control variables, while RCDs consist of finite panels so that the achievable magnitudes of solar propulsion and torque are actually discrete rather than continuous. This point was noticed by Aliasi et al. (2013) in their exploration of the active stabilization of L1-type artificial equilibrium points (AEPs) for a solar sail using RCDs. They considered the discretization and saturation effect and successfully generated AEPs by switching on or off the RCDs, which are symmetrically distributed on the sail. This paper naturally extends the application of finite discrete RCDs to the coupled attitude–orbit control of a spinning rigid solar sail in a Venus rendezvous mission. The paper aims to determine the switch law of RCDs, i.e., which part and when the devices should be switched on
to follow the desired orbit and attitude. The coupled attitude–orbit dynamical equation is derived from the SRP force model of RCDs implemented in IKAROS. A type of gradient projection algorithm is successfully developed to efficiently calculate the optimal switch sequences of RCDs, which is verified in the application to a Venus rendezvous mission. This paper are outlined as follows. In Section 2, the model of coupled attitude–orbit dynamics of spinning solar sail using RCDs is derived. Section 3 rearranges the derived linear form into an underdetermined linear system (ULS) with a binary set constraint (BSC), converts the formulated multi-objective least squares estimation (MLSE) and l1–l2 hybrid optimization into a classical quadratic programming (QP), and develops a type of basic gradient projection (BGP) method. As numerical verifications, switch regulations of RCDs are programmed applying previously proposed methods to accomplish the Venus rendezvous mission in Section 4. Section 5 draws conclusions of this paper. 2. Coupled attitude–orbit dynamics model of spinning solar sail utilizing reflectivity control devices 2.1. Configuration and coordinates of solar sail As depicted in Fig. 1, the solar sail is assumed to be an ultra-thin (zero thickness) rigid sheet with arbitrary contour shape and a total mass ms. To maintain the stability and deployment, the sail is spinning along its normal axis at a constant angular rate X. On the surface of the sail, RCDs are fixed with rpi as the vector from the mass center of the sail to the ith unit. A series of right-handed orthogonal reference coordinates are established. The J2000 ecliptic coordinates are chosen as the inertial coordinates oIxIyIzI to describe the
Fig. 1. Configuration (left) and coordinates (right) of solar sail.
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
orbital motion and the SRP force exerted on the sail. The sunlight pointing coordinates oExEyEzE are used to describe the relative attitude of the sail with respect to the sunlight. The origin oE is located at the sail’s mass center, the zE axis points from the Sun’s center to the sail’s center of mass, the yE axis is parallel to the orbital normal of the sail and the xE axis is located in the orbital plane and forms a right-handed coordinate system with yE and zE axes. Both the spin-free coordinates oFxFyFzF and the body-fixed coordinates oBxByBzB are established to describe the attitude motion and the SRP torque of the sail. The origins oF and oB are located at the same point as oE. The zF and zB axes are directed along the normal of the sail surface away from the Sun. Both the xF, xB and yF, yB axes are located in the sail plane; meanwhile, the xB axis rotates with the sail, and the xF axis does not follow the spinning motion. Assume that at the very beginning of the motion, the coordinates oFxFyFzF and oBxByBzB coincide with each other. Define the sun angle between the zF axis and the zE axis as a. The angle between the projection of the zF axis in the xEoEyE plane and the xE axis is denoted as d. As the sail is spinning along the normal vector of the sheet surface at constant angular speed X, the spin angle n can be easily obtained as n ¼ Xt
ð1Þ
where t denotes the spinning time. The intersection of the ecliptic plane and the orbital plane is defined as the node line oINI, and the angles between xI and oINI, and between oINI and zE are denoted as w and u, respectively. In addition, the angle h is the one between the yE and zI axes. The transition from oIxIyIzI to oExEyEzE and from oExEyEzE to oBxByBzB follows the 3-1-2 order with the angle sequence (w, p/2 + h, p/2 + u) and the 3-2-3 order with the angle sequence (d, a, n), respectively. Denote cosa and sina as ca and sa, respectively. The transition matrices from oIxIyIzI to oExEyEzE and from oExEyEzE to oBxByBzB can be represented as RIE ¼ Rz ðwÞRx ðp=2 þ hÞRy ðp=2 þ uÞ 2 3 cwsu swchcu swsh cwcu swchsu 6 7 ¼ 4 cwchcu swsu cwsh swcu þ cwchsu 5 shcu ch shsu ð2Þ and REB ¼ Rz ðdÞRy ðaÞRz ðnÞ 2 cacdcn sdsn cnsd cacdsn 6 ¼ 4 cdsn þ cacnsd cdcn casdsn cnsa
sasn
3 cdsa 7 sasd 5;
ð3Þ
ca
respectively. The transition matrix from oIxIyIzI to oBxByBzB can be derived as RIB = REBRIE. Angles w, h and u describe the orbital motion of the sail, and angles d, a and n represent the attitude motion of the sail.
1149
2.2. Force and torque models utilizing RCDs The surface of the sail is covered with k devices (k 2 N*) of RCDs in total. Each device can be switched between different reflection modes by powering on and off. Number the k devices from yB axis clockwise as i = 1, 2, . . ., k, with Dsi as the surface area of each device, and the SRP force Fi (McInnes, 1999; Rios-Reyes and Scheeres, 2005) exerted on the ith device can be written as 2 qd þ 2qs ðnR nÞ n F i ¼ P Dsi ðnR nÞ ðqa þ qd ÞnR þ 3 ð4Þ in which P represents the SRP on the sail, which can be calculated (Rios-Reyes and Scheeres, 2005) as 3 2pI 0 2 P¼ 1 ½1 ðRs =rÞ 2 ð5Þ 3c where I0 represents the frequency integrated specific intensity, c represents the speed of light, Rs represents the sun’s radius and r represents the distance from the sun; nR is the unit sun vector parallel to the zE axis; n is the unit normal vector of the sail parallel to the zF axis, and qa, qd and qs represent the coefficients of absorption, diffuse and specular reflections, respectively. In order to simplify the theoretical analysis, coefficients take values (Mu et al., 2015a,b) as qa_on = 0, qd_on = 0, qs_on = 1, qa_off = 0, qd_off = 1 and qs_off = 0, which is similar to the SRP model implemented in the IKAROS mission (Tsuda et al., 2011a, b). Thus, Eq. (4) can be rewritten as 2 F i ¼ P Dsi cos a nR þ n qdi þ 2 cos anqsi ð6Þ 3 in which coefficients qdiqsi = 0 and qdi + qsi = 1. If qdi is substituted by 1 qsi, qsi takes a value from the following binary set qsi 2 f0; 1g
ð7Þ
Thus, Eq. (6) can be rewritten as 2 F i ¼ P Dsi cos a 2 cos a n nR qsi 3 2 þ P Dsi cos a nR þ n 3
ð8Þ
Denote vector qs = [qs1, qs2, . . ., qsk]T and parameter Ps = P cosa. The total SRP force FR exerted on the sail can be obtained as k 2 2 X Dsi F R ¼ P s 2 cos a n nR S T qs þ P s nR þ n 3 3 i¼1 ð9Þ where qsi 2 {0, 1} and the vector S 2 R1 k is defined as S ¼ ½ Ds1
Ds2
Dsk
T
Define the matrix CF and the vector bF as
ð10Þ
1150
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
2 2 cos a n nR S T 3
respectively. The total SRP force can be rearranged as the following linear form
where r is the position vector from the Sun’s center oI to the center of mass of the solar sail oE and as is the SRP acceleration exerted on the sail as the dominant term of nonspherical geopotential J2 of the Sun is neglected. Since the performance of the J2 perturbation is much weaker (approximately 3 times less in order of magnitude) than SRP (Roxburgh, 2001), the neglect of J2 perturbation is reasonable. Derived from Eqs. (13) and (17), the orbit dynamics can be expressed in oIxIyIzI as
F R ¼ CF qs þ bF
ms as ¼ ms€r þ ðlms =r3 Þr ¼ F R þ eF ¼ I CF qs þ I bF þ eF
CF ¼ P s
ð11Þ
and k 2 X Dsi ; bF ¼ P s nR þ n 3 i¼1
ð12Þ
ð13Þ
ð18Þ
The torque Ti due to the SRP acting on the ith device can be obtained as 2 T i ¼ rpi F i ¼ P s Dsi rpi 2 cos a n nR qsi 3 2 ð14Þ þrpi nR þ n 3
where the superscript to the left of the symbol v means that the matrix or vector v is expressed in oIxIyIzI, and eF denotes the error between the expected solar force and the practical solar force generated by RCDs. The values of ICF and IbF can be calculated as
and the total torque of RCDs TR can be obtained as " # k X 2 Dsi rpi qsi T R ¼ nR 2 cos a n P s 3 i¼1 ! k X 2 Dsi rpi nR þ n P s : ð15Þ þ 3 i¼1
and
Denote the skew symmetric matrix of a vector v as v, and Eq. (15) can be rewritten as the following linear form 2 ½Dsi rpi 3k qs T R ¼ P s nR 2 cos a n 3 ! k X 2 þ Dsi rpi nR þ n P s ¼ CT qs þ bT ð16Þ 3 i¼1 in which the matrix [Dsirpi]3k = [Ds1rp1, Ds2rp2, . . ., Dskrpk]. 2.3. Coupled attitude–orbit dynamics The orbit dynamics can be expressed as the following perturbed equation in oIxIyIzI l €r ¼ 3 r þ as ð17Þ r
B
xEI
I
I
nR ¼ r=krk2
I
n ¼ REI RBE ½ 0
ð19Þ T
0 1
ð20Þ
where matrices RIE and REB can be obtained as Eqs. (2) and (3), respectively. Eq. (20) indicates that the orbital motion of the sail is coupled with the attitude motion due to the time-varying direction of vector n in the inertial coordinate. The angular velocity of the sail xBI can be obtained as xBI ¼ xBE þ xEI
ð21Þ
in which the angular velocity of oBxByBzB to oExEyEzE and oExEyEzE to oIxIyIzI can be expressed, respectively, in oBxByBzB as 2 3 2 3 2 3 0 0 0 6 7 6 7 6 7 T T T B xBE ¼ Rz ðnÞRy ðaÞ4 0 5 þ Rz ðnÞ4 a_ 5 þ 4 0 5 d_ 0 n_ 2 3 ð22Þ _ _ dsacn asn 6 7 _ ¼ 4 acn 5 _ þ dsasn _dca þ n_ and
2 3 2 39 0 > h_ = 6 7 6 7 6 7 T T T ¼ REB Ry ðp=2 þ uÞRx ðp=2 þ hÞ4 0 5 þ Ry ðp=2 þ uÞ4 0 5 þ 4 u_ 5 > > ; : 0 w_ 0 2 3 _ wcushÞ _ _ _ þ wsushÞðcdcn _ _ ðsdsn cacdcnÞðhsu ðcnsd þ cacdsnÞðu_ þ wchÞ þ cdsaðhcu casdsnÞðu_ þ wchÞ 6 7 _ wcushÞ _ _ þ wsushÞ _ ¼4 5 ðcdsn þ cacnsdÞðhsu þ sdsaðhcu 8 > <
2
0
3
_ þ wsushÞ _ _ _ wcushÞ _ caðhcu þ sasnðu_ þ wchÞ þ cnsaðhsu
ð23Þ
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
Matrices Rx, Ry and Rz are the rotating matrices along the x-, y- and z-axes, respectively. The superscript to the left of the symbol Bx means that the vector x is expressed in oBxByBzB. Denote the angular momentum of the sail as H = JBxBI, where the inertial moment matrix J = diag(Jx, Jy, Jz) is expressed in oBxByBzB and Jx, Jy and Jz and are the inertia moments about the xB, yB and zB axes, respectively. The attitude dynamics can be expressed using the following equation of the theorem of angular momentum B
_ þ B xBI H ¼ B T R þ eT TE ¼ H ¼ B CT qs þ B bT þ eT
s:t:qsi 2 f0; 1g
ð24Þ
where BTE represents the expected external torque expressed in oBxByBzB, and eT denotes the error between BTE and the torque generated from the RCDs. By substituting Eqs. (22) and (23) into Eq. (24), the equation can be rewritten in the form of Euler angles d, a and n. The unit vector along the yE axis and oINI expressed in oIxIyIzI is denoted as IyE and InI, respectively, which can be obtained as I
yE ¼ I r I v=kI r I vk2
ð25Þ
and I
nI ¼ ½ 0
0
ð26Þ
T
1 I yE
where vector Iv is the velocity of the sail expressed in oIxIyIzI. Furthermore, the orbital angles w, h and u can be obtained as 1 Þ
ð27Þ
0 0 T I nI Þ
ð28Þ
h ¼ arccosðI yE ½ 0 w ¼ arccosð½ 1
0
T
and u ¼ arccosðI nR I nI Þ
ð29Þ
By substituting BnR = REBRIEInR and Bn = REBRIEIn into Eq. (16), the solar torque BTR can be obtained and the theorem of angular momentum can finally be expressed in oBxByBzB. 3. Optimization for switch regulation of RCDs 3.1. Underdetermined Linear System with Binary Set Constraint The switch programming aims to attain an appropriate switch sequence for RCDs to achieve the desired control
qs
"
force and torque. Denote the programming vector Y = [msas, BTE]T 2 R61, and the equation of coupled attitude–orbit dynamics can be written as the following linear form Y ¼ Cqs þ b þ e
s:t:qsi 2 f0; 1g
ð30Þ
where matrix C, vector b and vector e can be represented as "I # "I # CF bF 6k C¼ 2R ; b¼ 2 R61 ; B B CT bT ð31Þ " # eF e¼ 2 R61 eT Naturally, the optimal value q*s of the power vector is sought to minimize the orbit and attitude errors, which can be expressed as the following multi-objective LSE (MLSE) problem 1 2 2 qs ¼ argqsi 2f0;1g min ðkeF k2 þ kT keT k2 Þ 2
ð32Þ
in which the denotation || ||2 represents the 2-norm of the vector, and kT is the Lagrange multiplier for weighting the attitude error. However, taking the electric energy consumption into account, it is expected to utilize as few switched-on RMDs as possible, if only the control accuracy is guaranteed; as a result, an l1–l2 hybrid optimization is developed as 1 kT 2 2 keF k2 þ keT k2 þ k1 kqs k1 qs ¼ argqsi 2f0;1g min ð33Þ 2 2 where the denotation || ||1 represents the 1-norm of the vector, and k1 is the Lagrange multiplier weighting the electric energy consumption. Furthermore, to provide the MLSE and the l1–l2 hybrid optimization with a unified quadratic programming (QP) scheme shown as follows, 1 T z Hz þ cT z z ¼ arg min zi 2f0;1g;16i6k 2 zj 2f1;0g;kþ16j62k
¼ arg
min
zi 2f0;1g;16i6k zj 2f1;0g;kþ16j62k
GðzÞ
ð34Þ
new matrices and vectors are defined as
ðCF þ kT CT ÞT ðCF þ kT CT Þ ðCF þ kT CT ÞT ðCF þ kT CT Þ
#
; ;H¼ T T qs ðCF þ kT CT Þ ðCF þ kT CT ÞðCF þ kT CT Þ ðCF þ kT CT Þ # " 1k1 CTF PF kT CTT PT c ¼ k1 þ 1k1 CTF PF þ kT CTT PT z¼
1151
ð35Þ
1152
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
Note that every element of power vector qs takes a value from the binary set {0, 1}, which is called the BSC. In addition, in most situations (when k P 6) the matrix C is underdetermined and does not have a Moore–Penrose inverse (a kind of pseudo inverse matrix, see (Zhang, 2004) for details). Therefore, the switch law design of RCDs is actually an underdetermined QP with BSC and has no analytical solution. As a result, a numerical method for the optimal search is sought. The problem has been converted into a QP in the previous discussion. Nevertheless, the BSC leads to an obstacle for directly applying the widely available algorithms (such as interior-point method, trust-region-reflectivity method, active-set method and so on) which rely on the first-order essential condition (i.e. 5G(z*) = 0) of the optimal solution, so the relaxation of the BSC to a continuous interval is performed before the optimal search. Meanwhile, a complete relaxation of BSC to the real number field will cause numerical solutions converging to an improper one. Thus, the BSC should be relaxed to the closed interval [1, +1] and an optimization problem can be derived from Eq. (34) as zs ¼ arg
min
zi 2½0;1;16i6k zj 2½1;0;kþ16j62k
GðzÞ
ð36Þ
in which the vector zs is a sub-optimal solution of problem (34). Furthermore, an approximate solution of problem (34) can be obtained after the discretization of zsi to the binary set {0, 1}, which can be written as ( 0 zsi < 0:5 qs ¼ ½qsi k1 ; qsi ¼ ; i ¼ 1; . . . ; k : 1 zsi P 0:5 ð37Þ
Fig. 3. The procedure of the BGP method.
3.2. Basic gradient projection: the BGP method A type of basic gradient projection (BGP) method (Bertsekas, 1999; Ma´rio et al., 2007) is proposed for the numerical calculation of Eq. (36). The geometry illustration of the BGP in a 2-dimensional case is given by Fig. 2. The BGP achieves the optimal search along the gradient descend of function G(z) with the complete relaxation of the BSC to the region of real number. After every step of descend, the sub-optimal solution is projected to the finite interval [1, +1] and the continuous optimal solution zs can be obtained. Subsequently, make the discretization of zsi to the binary set {0, 1} and finally obtain the optimal solution q*s . Denote the projection operator Proj as Proj½v ¼ ½vf ¼ ½vi k1 ;
Fig. 2. The geometry illustration of the BGP in a 2-dimensional case.
vi ¼ midð0; vi ; 1Þ
ð38Þ
in which function mid(a, b, c) denotes the middle value of its three scalar arguments. Define r as the damping factor of the backtracking, h as the backtracking step of argument, and l as the backtracking step of function value. The procedure of the BGP method is depicted in Fig. 3 and can be presented as follows. k1 Step 0 (initialization): Choose q(0) | s 2 W = {qs 2 R qsi 2 [1, +1]}; r 2 (0, 1), l 2 (0, 1/2) and [hmin, hmax] satisfying 0 < hmin. For p = 0, 1, 2, . . . Step 1: Compute h(0): h(0) = ((g(p))Tg(p))/((Hg(p))THg(p)). Replace h(0) by mid(hmin, h(0), hmax).
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
1153
Table 1 Parameters of the sail. Parameter
Specification
Mass (ms) Sail radius (rs) Payload region radius (rl) Interior SRP radius (rpin) Exterior SRP radius (rpex) Interior unit surface area (DSin) Exterior unit surface area (DSex) Moments of inertia (Jx, Jy, Jz) Spin rate (X)
64 kg 30 m 1m 9.3 m 29.5 m 77 m2 1.5 m2 (14,400, 14,400, 28,800) kg m2 2.5 rpm
Table 2 Parameters of the Venus rendezvous mission.
Fig. 4. The structure of solar sail demo.
Step 2 (backtracking line search): Choose h(p) = h(0), (0) rh , r2h(0), . . . such that (p) (p) (p) (p) T (p) G([q(p) s h 5G(qs )]f) 6 G(qs ) l[5G(qs )] {qs (p) (p+1) (p) (p) [q(p) = [q(p) 5G s h 5G(qs )]f} and set qs s h (q(p) s )]f. Step 3 (termination test): End the search if the termination condition is satisfied; otherwise, return to Step 1 for a new iteration. Here, the termination condition is chosen as kGðqsðpþ1Þ Þ GðqðpÞ s Þk2 6 tolG
ð39Þ
Vector g is defined as ( ðpÞ ðrGðqðpÞ qsi > 0 or ðrGðqðpÞ ðpÞ s ÞÞi ; s ÞÞi < 0 gi ¼ 0; otherwise: (p)
ð40Þ 4. Numerical demonstration The coupled attitude–orbit dynamics of a solar sail utilizing RCDs is numerically calculated in a Venus rendezvous mission. This section will present and analyse numerical results of both MLSE and l1–l2 hybrid optimization given by Eqs. (32) and (33). 4.1. Sail and mission description A circular solar sail depicted as Fig. 4 is used as a demo for detailed numerical calculations. The sail is a rigid even circular sheet with radius rs and total mass ms (payload mass included). On the surface of the sail, there are 16 independently switched RCDs in which 8 units are fixed on the interior region where rpin is the radius from oB to the interior SRP center and DSin is the surface area for every interior device; meanwhile, the other 8 units are fixed on the
Parameter
Specification
Initial time (TI) Initial position (xI, yI, zI) Initial velocity (vxI, vyI, vzI) Terminal time (TT) Terminal position (xT, yT, zT) Terminal velocity (vxT, vyT, vzT) Transfer time (tT)
24 May 2016 04:19:03.360 UTCG (0.4622, 0.9032, 1.4647 106) AU (2.6008 104, 1.3591 104, 2.2335) m/s 12 Sep 2017 05:07:17.76 UTCG (0.0874, 0.7196, 0.0149) AU (3.4599 104, 4.4822 103, 1.9295 103) m/s 476.0335 day
edge of the sail where rpex is the radius from oB to the exterior SRP center and DSex is the surface area for every interior device. Surrounded by the interior RCD, a circular region with radius rl is designed to have payloads (e.g., satellite computer, camera, etc.) fastened on. The remaining region of the sail is mainly membrane to hold the structure and collect solar energy. Table 1 presents detailed values for parameters of the sail. Notice that DSi is much larger than DSe, while rpi is much smaller than rpe. Thus, the solar acceleration is mainly generated by interior RCDs, while the torque is mainly generated by the exterior RCDs. The attitude control and the orbit control can be partly decoupled in structure. Besides, the circular sheet configuration is beneficial for the on-orbital deployment and the location of payloads will result in small centrifugal force or turbulent torque. For the Venus rendezvous, it is assumed that the solar sail departs from the Earth with no escape velocity, i.e., C3 = 0 and no propulsion occurs at the end of the mission. The reference orbit and attitude are calculated using an indirect method (Gong and Li, 2014; Jayaraman, 1980; Wood, 1982) to optimize the transfer time, assuming that the reflectivity modulation ratio is continuous. Parameters of the Venus rendezvous mission are presented in Table 2. 4.2. Numerical results and analysis Consider the MLSE problem given by Eq. (33) to ensure the accuracy of orbit control as well as compensate some
1154
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
Terminal Position Earth zI ˄AU˅
Reference Orbit Practical Orbit
Initial Position
Venus
xI˄AU˅ yI˄AU˅
Fig. 5. Practical transfer orbit of the sail when kT = 0.5.
5
2
x 10
100
1
0
50
0
50
100
150
200
250 300 Time(day)
350
400
450
Velocity Error(m/s)
Position Error(km)
Position Error Velocity Error
0 500
Fig. 6. Position and velocity error of the transfer orbit when kT = 0.5.
part of desired torque. Set the Lagrange multiplier kT = 0.5 for weighting attitude error. The numerical results are illustrated as Figs. 5–10. Fig. 5 illustrates the reference and practical transfer orbits which coincide with each other in the figure and can hardly be distinguished. Therefore, it is necessary to evaluate the ability of the switch programming by Fig. 6 which illustrates increasing errors in position and velocity that finally reach approximately 105 km and 50 m/s, respectively. The order of the terminal position error indicates that the sail finally locates beneath the influence surface of the Venus; meanwhile the terminal velocity of the sail is close to the Venus. Thus, the numerical result is compatible with the reference mission. These errors can be reduced and finally eliminated by increasing the total amount and distributed surface area of RCDs and decreasing the surface area of every single RCD, as indicated by Fig. 14 and Mu et al. (2013, 2015a,b). Fig. 7 shows the expected solar acceleration and practical one generating from RCDs along three axes of the inertial coordinate. The SRP
Fig. 7. SRP acceleration and error when kT = 0.5.
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
1155
Fig. 8. SRP torque compensation when kT = 0.5.
acceleration is approximately 104 m/s2, while the error between the expected and practical SRP acceleration is approximately 107 m/s2. Both of the expected and practical SRP torques along three axes of the body-fixed coordinate are illustrated in Fig. 8, which indicates a partial compensation of the torque by utilizing RCDs. However, as the SRP force given by Eq. (4) has a component along vector nR, the component of the SRP torque along zB axis is approximately 103 Nm, which performs as a disturbance of the spinning stability. Figs. 9 and 10 illustrate the discrete switch law for the eight interior RCDs and the eight exterior RCDs, respectively. Clearly, the interior RCDs mainly contribute to the orbit control and are switched on for nearly the entire voyage; while the exterior RCDs mainly contribute to the attitude control and are switched between the on and off states much more
frequently. Furthermore, as the Lagrange multiplier kT varies, the optimization objective and numerical consequences of switch programming change, which are illustrated in Figs. 11–13. Define the maximum relative error of SRP acceleration as ea max ¼ maxðkeF k2 =kms as k2 Þ 100%
ð41Þ
where function max(s) represents the maximum value of the scalar s. Illustrated by Figs. 11 and 12, as kT increases from 0 to 1, eamax increases from approximately 0.03% to 0.055%; while terminal position and velocity error increases from approximately 1.5 105 km to 1.65 105 km and from approximately 48 m/s to 53 m/s, respectively. On the other hand, define the uncompensated time of SRP torque as
Fig. 9. Discrete switch law for interior RCDs when kT = 0.5.
1156
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158 1
Unit 1 0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0
50
100
150
200
250
300
350
400
450
500
Unit 2 Unit 3 Unit 4 Unit 5 Unit 6 Unit 7 Unit 8
Time (day)
Fig. 10. Discrete switch law for exterior RCDs when kT = 0.5.
0.055
0.05
400 350
0.045
Uncompensated time(day)
Maximum relative error of acceleration(%)
450
0.04
0.035
0.03
300 250 200 150 100
0.025
0
0.1
0.2
0.3
0.4
0.5 0.6 Lambda
0.7
0.8
0.9
1
50
0
0.1
0.2
0.3
0.4
Fig. 11. Maximum relative error of SRP acceleration versus kT.
0.5 0.6 Lambda
0.7
0.8
0.9
1
Fig. 13. Uncompensated time of SRP torque utilizing RCD versus kT.
5.5
x 10
Terminal Position Error Terminal Velocity Error
TerminalPosition Error Terminal Velocity Error
1.6
52
1.5
50
Position Error(log(km))
5
Velocity Error(m/s)
Position Error(km)
50
54
45 40 35
4.5 30 4
25 20
3.5
Velocity Error(m/s)
5
1.7
15 10
3 5
1.4 0
0.1
0.2
0.3
0.4
0.5 0.6 Lambda
0.7
0.8
0.9
48 1
Fig. 12. Terminal position and velocity error of rendezvous versus kT.
2.5 10
20
30
40
50
60
70
0 80
k
Fig. 14. Terminal position and velocity error of the Venus rendezvous versus k.
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
1157 80
1500
80
L1=0,Position Error L1=0,Switched-on RMD
1000
Postion Error(km) 500
79 Switched-on 78 RMD
79
78
00 00
50 50
100 100
150 150
250 250
200 200 200
300 300 300
350 350 350
400 400 400
450 450
10000
77 77 500 500 500 80
80
Postion Error(km) 5000
Switched-on 78 RMD
78
L1=5e-6,Position Error L1=5e-6,Switched-on RMD 00 00
50 50 50
100 100 100
150 150 150
250 250 250
200
200 200
300 300 300
350
350 350
400 400 400
450
450 450
76 76 500 500 500
5
2
Postion Error(km)
x 10
100
100
L1=1e-5,Position Error L1=1e-5,Switched-on RMD
Switched-on 50 RMD
1
00 00
50
50 50
100 100
150 150
200 200 200
250 250 Time (day)
300 300 300
350 350 350
400 400 400
450 450
0 0 500 500 500
Fig. 15. Position error and number of switched-on RCD unit for l1–l2 hybrid optimization.
tunc ¼ timeðkeT k2 =kB T E k2 > eT Þ
ð42Þ
where function time(S) represents the total time that the statement S is true, and eT represents the threshold of the truth of statement S. Here, eT is set to equal 0.5. It is illustrated in Fig. 13 that the uncompensated time of SRP torque decreases from approximately 400 days to 100 days when kT increases from 0 to 1. Numerical results in Figs. 11–13 clearly demonstrate the weighting effect of the Lagrange multiplier for MLSE. As the value of kT increases, the optimal switch law of RCDs strengthens the compensation of desired torque and weakens the accuracy of orbit control, which can also be deduced from the mathematical expression of MLSE shown as Eq. (32). Since the reference orbit and attitude are obtained by assuming that the modulation surface area of RCDs is continuous, a decrease of both terminal position and velocity errors is expected with an increase in the total amount of RCDs. This inference can be verified numerically in the Venus rendezvous mission. Assuming kT = 0 which means the MLSE only considers the orbital control accuracy, recycle the switch programming of RCDs for different k values (the total amount of RCDs) and calculate the terminal position and velocity errors. The numerical result is shown in Fig. 14 in which the left vertical axis ‘Position Error (log(km))’ represents the log10 value of the terminal position error. It is clear that versus the increase of k from 16 to 80, the terminal position error falls from approximately 10 km to 103 km and the terminal velocity error falls from approximately 50 m/s to no more than 1 m/s. Thus, as stated previously, it is reasonable to infer that the orbital control error will finally be eliminated by RCDs as the modulation ratio tends towards a continuous variable. Finally, the l1–l2 hybrid optimization defined as Eq. (33) is computed numerically. Set k = 80 and kT = 0.01 and consider three situations: k1 = 0, k1 = 5 106 and k1 = 105. The position error and the number of switched-on
RCDs are illustrated in Fig. 15. The increase of the Lagrange multiplier k1 adds the weight of the consumption of electrical energy for the optimization. This effect is clearly illustrated in the numerical results as an obvious increase of terminal position error from approximately 1500 km to 105 km as well as a decrease of amount of switched-on RCDs. 5. Conclusions This paper derived coupled attitude–orbit dynamics of a solar sail utilizing independently switched reflectivity control devices (RCDs). According to the derivation, the solar force and torque were linearly dependent on the switch vector qs which usually had a higher dimension than the summing dimension of the force and torque. In addition, every element of vector qs took a value from a binary set {0, 1} due to the switch property of RCDs. Thus, the coupled attitude–orbit dynamics could be described by an underdetermined linear system (ULS) subject to the binary set constraint (BSC) in mathematics. This discovery helped the authors formulate the switch programming of RCDs into a multi-objective least-squares estimation (MLSE) and an l1–l2 hybrid optimization, respectively. By reformulating these two optimization problems into a unified quadratic programming (QP), a basic gradient projection (BGP) method was designed for the search of optimal solutions. A circular rigid sail with RMDs flying in the Venus rendezvous mission was taken as a demonstration of the derived dynamics and proposed BGP. Numerical results indicated that with an appropriate switch law of RMDs, the sail achieved to arrive beneath the influence surface of the Venus with approximately 105 km position error and 50 m/s velocity error. RMDs also generated solar torque which could compensate some part of desired angular acceleration of the sail. Furthermore, numerical results also illustrated that the terminal position and velocity error decreased below approximately 103 km and 0.5 m/s,
1158
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
respectively, as the total amount of RCDs increased to 80. This phenomenon indicated a promising way to reduce errors for orbital control of the sail. The switch programming of RCDs in this paper was a kind of open-loop control, so a feedback controller was expected to eliminate the errors in the future. On the other hand, every element of the switch vector took a value from the binary set {0, 1}, which was a simplification for theoretical analysis. As a matter of fact, values of optical coefficients varied as the wavelength, and accurate synthetic values was looked forward to by on-orbital experiments. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grants no. 11272004). The authors wish to thank Fangming Han from Tsinghua University for his introduction to the underdetermined linear system with a binary set constraint. They also appreciate Pascal Willis and both anonymous reviewers for their comments and kind help. References Aliasi, G., Mengali, G., Quarta, A., 2013. Artificial Lagrange points for solar sail with electrochromic material panels. J. Guidance Control Dyn. 36, 1544–1550. Bertsekas, D., 1999. Nonlinear Programming, Athena. Gong, S., Baoyin, H., Li, J., 2009. Coupled attitude–orbit dynamics and control for displaced solar orbits. Acta Astronaut. 65, 730–737. Gong, S., Li, J., 2014. Fuel consumption for interplanetary missions of solar sailing. Sci. China-Phys. Mech. Astron. 57, 521–531. Jayaraman, T., 1980. Time-optimal orbit transfer trajectory for solar sail spacecraft. J. Guidance Control Dyn. 3, 536–542. Macdonald, M., McInnes, C., 2011. Solar sail science mission applications and advancement. Adv. Space Res. 48, 1702–1716. Ma´rio, A., Robert, D., Stephen, J., 2007. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems, IEEE. J. Sel. Top. Signal Process. 1, 586–597.
McInnes, C.R., 1999. Solar Sailing: Technology, Dynamics and Mission Applications. Springer-Verlag, Berlin (pp. 90–95). Mu, J., Gong, S., Li, J., 2013. Reflectivity-controlled solar sail formation flying for magneto sphere mission. Aerosp. Sci. Technol. 30, 339–348. Mu, J., Gong, S., Li, J., 2015a. Coupled control of reflectivity modulated solar sail for GeoSail formation flying. J. Guidance Control Dyn. 38, 740–751. Mu, J., Gong, S., Ma, P., Li, J., 2015b. Dynamics and control of flexible spinning solar sails under reflectivity modulation. Adv. Space Res. 56, 1731–1751. Rios-Reyes, L., Scheeres, D., 2005. Generalized model for solar sails. J. Spacecraft Rockets 42, 182–185. Roxburgh, I., 2001. Gravitational multipole moments of the Sun determined from helioseismic estimates of the internal structure and rotation. A. & A. 377, 688–690. Scholz, C., Romagnoli, D., Dachwald, B., Theil, S., 2011. Performance analysis of an attitude control system for solar sails using sliding masses. Adv. Space Res. 48, 1822–1835. Tsuda, Y., Mori, O., Funase, R., Sawada, H., et al., 2011a. Flight status of IKAROS deep space solar sail demonstrator. Acta Astronaut. 69, 833– 840. Tsuda, Y., Saiki, T., Funase, R., Mimasu, Y., 2011b. On-orbit verification of fuel-free attitude control system for spinning solar sail utilizing solar radiation pressure. Adv. Space Res. 48, 1740–1746. Tsuda, Y., Saiki, T., Funase, R., Mimasu, Y., 2013. Generalized attitude model for spinning solar sail spacecraft. J. Guidance Control Dyn. 36, 967–974. Wie, B., 2004a. Solar sail attitude control and dynamics, Part 1. J. Guidance Control Dyn. 27, 526–535. Wie, B., 2004b. Solar sail attitude control and dynamics, Part 2. J. Guidance Control Dyn. 27, 536–544. Wood, L., Bauer, T., Zondervan, K., 1982. Comment on ‘‘time-optimal orbit transfer trajectory for solar sail spacecraft”. J. Guidance Control Dyn. 5, 221–224. Zeng, X., Alfriend, K.T., Vadali, S.R., 2014. Solar sail planar multireversal periodic orbits. J. Guidance Control Dyn. 37, 674–681. Zhang, J., Wang, T., 2013. Coupled attitude–orbit control of a flexible solar sail for a displaced solar orbit. J. Spacecraft Rockets 50, 675–685. Zhang, X., 2004. Matrix Analysis and Applications. Tsinghua University, Beijing (pp. 85–96).
ScienceDirect Advances in Space Research 57 (2016) 1147–1158 www.elsevier.com/locate/asr
Switch programming of reflectivity control devices for the coupled dynamics of a solar sail Tianjian Hu a,b,⇑, Shengping Gong a, Junshan Mu a,c, Junfeng Li a, Tianshu Wang a, Weiping Qian b b
a School of Aerospace Engineering, Tsinghua University, Beijing 100084, China Beijing Institute of Tracking and Telecommunication Technology, Beijing 100094, China c China Satellite Maritime Tracking and Control Department, Jiangyin 214431, China
Received 15 September 2015; received in revised form 16 December 2015; accepted 17 December 2015 Available online 25 December 2015
Abstract As demonstrated in the Interplanetary Kite-craft Accelerated by Radiation Of the Sun (IKAROS), reflectivity control devices (RCDs) are switched on or off independently with each other, which has nevertheless been ignored by many previous researches. This paper emphasizes the discrete property of RCDs, and aims to obtain an appropriate switch law of RCDs for a rigid spinning solar sail. First, the coupled attitude–orbit dynamics is derived from the basic solar force and torque model into an underdetermined linear system with a binary set constraint. Subsequently, the coupled dynamics is reformulated into a constrained quadratic programming and a basic gradient projection method is designed to search for the optimal solution. Finally, a circular sail flying in the Venus rendezvous mission demonstrates the model and method numerically, which illustrates approximately 103 km terminal position error and 0.5 m/s terminal velocity error as 80 independent RCDs are switched on or off appropriately. Ó 2015 COSPAR. Published by Elsevier Ltd. All rights reserved.
Keywords: Solar sail; Coupled attitude–orbit dynamics; Reflectivity control devices; Gradient projection method
1. Introduction Due to its capability of generating propulsion without fuel consumption, the solar sail has been envisioned to enable or enhance a wide range of mission applications throughout the solar system (Macdonald and McInnes, 2011). Some previous literature (Gong and Li, 2014; Zeng et al., 2014) mainly focused on the orbit design of a solar sail and assumed that the attitude can change
⇑ Corresponding author at: School of Aerospace Engineering, Tsinghua University, Beijing 100084, China. E-mail addresses: [email protected] (T. Hu), [email protected] (S. Gong), [email protected] (J. Mu), [email protected]. cn (J. Li), [email protected] (T. Wang), qianweipingbittt@sina. com (W. Qian).
http://dx.doi.org/10.1016/j.asr.2015.12.029 0273-1177/Ó 2015 COSPAR. Published by Elsevier Ltd. All rights reserved.
instantaneously. However, the propulsion exerted by solar radiation pressure (SRP) is related to the angle between the sunlight and the normal vector of the sail surface, so the orbit dynamics and attitude dynamics are not completely independent of each other (Gong et al. 2009; Zhang and Wang, 2013). Thus, neither the attitude maneuver process nor the attitude–orbit coupled effects can be neglected. Starting from the above idea, Gong et al. (2009) developed the coupled attitude–orbit dynamics for displaced solar orbits and discussed the stability of its motion. Zhang and Wang (2013) considered the flexibility of the sail membrane and derived a coupled dynamical equation for a flexible displaced solar orbiter. Nevertheless, neither of these authors discussed the physical realization of the desired trajectories in both orbit and attitude. An interesting implementation of reflectivity control devices (RCDs)
1148
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
on the attitude adjustment of a solar sail was demonstrated by the Japanese sail ‘‘interplanetary kitecraft accelerated by radiation of the sun” (IKAROS) (Tsuda et al., 2011a,b, 2013). RCD is a thin film attached to the surface of the sail that can switch its reflectivity between the specular and diffusive states versus electricity power being switched on or off. Different reflectivity states can result in different SRP. By adjusting the total amount and distribution of switched-on RCDs, the propulsion and torque acting on the sail can be controlled to follow the desired orbit and attitude. Compared to other methods for attitude control utilizing the gimbaled control boom (Wie 2004a,b), sail control vanes (Wie 2004a,b), or sliding masses (Scholz et al., 2011), the application of RCDs is obviously propellantless, much simpler in structure and has lower cost (Tsuda et al., 2011a,b). Such benefits have prompted increasing research interest in the utilization of RCDs in solar sail missions. Mu et al. (2013, 2015a,b) extended the coupled attitude–orbit dynamics to formation flying for a GeoSail mission. The reflectivity modulation ratio and its distribution on the sail are used as continuous control variables, while RCDs consist of finite panels so that the achievable magnitudes of solar propulsion and torque are actually discrete rather than continuous. This point was noticed by Aliasi et al. (2013) in their exploration of the active stabilization of L1-type artificial equilibrium points (AEPs) for a solar sail using RCDs. They considered the discretization and saturation effect and successfully generated AEPs by switching on or off the RCDs, which are symmetrically distributed on the sail. This paper naturally extends the application of finite discrete RCDs to the coupled attitude–orbit control of a spinning rigid solar sail in a Venus rendezvous mission. The paper aims to determine the switch law of RCDs, i.e., which part and when the devices should be switched on
to follow the desired orbit and attitude. The coupled attitude–orbit dynamical equation is derived from the SRP force model of RCDs implemented in IKAROS. A type of gradient projection algorithm is successfully developed to efficiently calculate the optimal switch sequences of RCDs, which is verified in the application to a Venus rendezvous mission. This paper are outlined as follows. In Section 2, the model of coupled attitude–orbit dynamics of spinning solar sail using RCDs is derived. Section 3 rearranges the derived linear form into an underdetermined linear system (ULS) with a binary set constraint (BSC), converts the formulated multi-objective least squares estimation (MLSE) and l1–l2 hybrid optimization into a classical quadratic programming (QP), and develops a type of basic gradient projection (BGP) method. As numerical verifications, switch regulations of RCDs are programmed applying previously proposed methods to accomplish the Venus rendezvous mission in Section 4. Section 5 draws conclusions of this paper. 2. Coupled attitude–orbit dynamics model of spinning solar sail utilizing reflectivity control devices 2.1. Configuration and coordinates of solar sail As depicted in Fig. 1, the solar sail is assumed to be an ultra-thin (zero thickness) rigid sheet with arbitrary contour shape and a total mass ms. To maintain the stability and deployment, the sail is spinning along its normal axis at a constant angular rate X. On the surface of the sail, RCDs are fixed with rpi as the vector from the mass center of the sail to the ith unit. A series of right-handed orthogonal reference coordinates are established. The J2000 ecliptic coordinates are chosen as the inertial coordinates oIxIyIzI to describe the
Fig. 1. Configuration (left) and coordinates (right) of solar sail.
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
orbital motion and the SRP force exerted on the sail. The sunlight pointing coordinates oExEyEzE are used to describe the relative attitude of the sail with respect to the sunlight. The origin oE is located at the sail’s mass center, the zE axis points from the Sun’s center to the sail’s center of mass, the yE axis is parallel to the orbital normal of the sail and the xE axis is located in the orbital plane and forms a right-handed coordinate system with yE and zE axes. Both the spin-free coordinates oFxFyFzF and the body-fixed coordinates oBxByBzB are established to describe the attitude motion and the SRP torque of the sail. The origins oF and oB are located at the same point as oE. The zF and zB axes are directed along the normal of the sail surface away from the Sun. Both the xF, xB and yF, yB axes are located in the sail plane; meanwhile, the xB axis rotates with the sail, and the xF axis does not follow the spinning motion. Assume that at the very beginning of the motion, the coordinates oFxFyFzF and oBxByBzB coincide with each other. Define the sun angle between the zF axis and the zE axis as a. The angle between the projection of the zF axis in the xEoEyE plane and the xE axis is denoted as d. As the sail is spinning along the normal vector of the sheet surface at constant angular speed X, the spin angle n can be easily obtained as n ¼ Xt
ð1Þ
where t denotes the spinning time. The intersection of the ecliptic plane and the orbital plane is defined as the node line oINI, and the angles between xI and oINI, and between oINI and zE are denoted as w and u, respectively. In addition, the angle h is the one between the yE and zI axes. The transition from oIxIyIzI to oExEyEzE and from oExEyEzE to oBxByBzB follows the 3-1-2 order with the angle sequence (w, p/2 + h, p/2 + u) and the 3-2-3 order with the angle sequence (d, a, n), respectively. Denote cosa and sina as ca and sa, respectively. The transition matrices from oIxIyIzI to oExEyEzE and from oExEyEzE to oBxByBzB can be represented as RIE ¼ Rz ðwÞRx ðp=2 þ hÞRy ðp=2 þ uÞ 2 3 cwsu swchcu swsh cwcu swchsu 6 7 ¼ 4 cwchcu swsu cwsh swcu þ cwchsu 5 shcu ch shsu ð2Þ and REB ¼ Rz ðdÞRy ðaÞRz ðnÞ 2 cacdcn sdsn cnsd cacdsn 6 ¼ 4 cdsn þ cacnsd cdcn casdsn cnsa
sasn
3 cdsa 7 sasd 5;
ð3Þ
ca
respectively. The transition matrix from oIxIyIzI to oBxByBzB can be derived as RIB = REBRIE. Angles w, h and u describe the orbital motion of the sail, and angles d, a and n represent the attitude motion of the sail.
1149
2.2. Force and torque models utilizing RCDs The surface of the sail is covered with k devices (k 2 N*) of RCDs in total. Each device can be switched between different reflection modes by powering on and off. Number the k devices from yB axis clockwise as i = 1, 2, . . ., k, with Dsi as the surface area of each device, and the SRP force Fi (McInnes, 1999; Rios-Reyes and Scheeres, 2005) exerted on the ith device can be written as 2 qd þ 2qs ðnR nÞ n F i ¼ P Dsi ðnR nÞ ðqa þ qd ÞnR þ 3 ð4Þ in which P represents the SRP on the sail, which can be calculated (Rios-Reyes and Scheeres, 2005) as 3 2pI 0 2 P¼ 1 ½1 ðRs =rÞ 2 ð5Þ 3c where I0 represents the frequency integrated specific intensity, c represents the speed of light, Rs represents the sun’s radius and r represents the distance from the sun; nR is the unit sun vector parallel to the zE axis; n is the unit normal vector of the sail parallel to the zF axis, and qa, qd and qs represent the coefficients of absorption, diffuse and specular reflections, respectively. In order to simplify the theoretical analysis, coefficients take values (Mu et al., 2015a,b) as qa_on = 0, qd_on = 0, qs_on = 1, qa_off = 0, qd_off = 1 and qs_off = 0, which is similar to the SRP model implemented in the IKAROS mission (Tsuda et al., 2011a, b). Thus, Eq. (4) can be rewritten as 2 F i ¼ P Dsi cos a nR þ n qdi þ 2 cos anqsi ð6Þ 3 in which coefficients qdiqsi = 0 and qdi + qsi = 1. If qdi is substituted by 1 qsi, qsi takes a value from the following binary set qsi 2 f0; 1g
ð7Þ
Thus, Eq. (6) can be rewritten as 2 F i ¼ P Dsi cos a 2 cos a n nR qsi 3 2 þ P Dsi cos a nR þ n 3
ð8Þ
Denote vector qs = [qs1, qs2, . . ., qsk]T and parameter Ps = P cosa. The total SRP force FR exerted on the sail can be obtained as k 2 2 X Dsi F R ¼ P s 2 cos a n nR S T qs þ P s nR þ n 3 3 i¼1 ð9Þ where qsi 2 {0, 1} and the vector S 2 R1 k is defined as S ¼ ½ Ds1
Ds2
Dsk
T
Define the matrix CF and the vector bF as
ð10Þ
1150
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
2 2 cos a n nR S T 3
respectively. The total SRP force can be rearranged as the following linear form
where r is the position vector from the Sun’s center oI to the center of mass of the solar sail oE and as is the SRP acceleration exerted on the sail as the dominant term of nonspherical geopotential J2 of the Sun is neglected. Since the performance of the J2 perturbation is much weaker (approximately 3 times less in order of magnitude) than SRP (Roxburgh, 2001), the neglect of J2 perturbation is reasonable. Derived from Eqs. (13) and (17), the orbit dynamics can be expressed in oIxIyIzI as
F R ¼ CF qs þ bF
ms as ¼ ms€r þ ðlms =r3 Þr ¼ F R þ eF ¼ I CF qs þ I bF þ eF
CF ¼ P s
ð11Þ
and k 2 X Dsi ; bF ¼ P s nR þ n 3 i¼1
ð12Þ
ð13Þ
ð18Þ
The torque Ti due to the SRP acting on the ith device can be obtained as 2 T i ¼ rpi F i ¼ P s Dsi rpi 2 cos a n nR qsi 3 2 ð14Þ þrpi nR þ n 3
where the superscript to the left of the symbol v means that the matrix or vector v is expressed in oIxIyIzI, and eF denotes the error between the expected solar force and the practical solar force generated by RCDs. The values of ICF and IbF can be calculated as
and the total torque of RCDs TR can be obtained as " # k X 2 Dsi rpi qsi T R ¼ nR 2 cos a n P s 3 i¼1 ! k X 2 Dsi rpi nR þ n P s : ð15Þ þ 3 i¼1
and
Denote the skew symmetric matrix of a vector v as v, and Eq. (15) can be rewritten as the following linear form 2 ½Dsi rpi 3k qs T R ¼ P s nR 2 cos a n 3 ! k X 2 þ Dsi rpi nR þ n P s ¼ CT qs þ bT ð16Þ 3 i¼1 in which the matrix [Dsirpi]3k = [Ds1rp1, Ds2rp2, . . ., Dskrpk]. 2.3. Coupled attitude–orbit dynamics The orbit dynamics can be expressed as the following perturbed equation in oIxIyIzI l €r ¼ 3 r þ as ð17Þ r
B
xEI
I
I
nR ¼ r=krk2
I
n ¼ REI RBE ½ 0
ð19Þ T
0 1
ð20Þ
where matrices RIE and REB can be obtained as Eqs. (2) and (3), respectively. Eq. (20) indicates that the orbital motion of the sail is coupled with the attitude motion due to the time-varying direction of vector n in the inertial coordinate. The angular velocity of the sail xBI can be obtained as xBI ¼ xBE þ xEI
ð21Þ
in which the angular velocity of oBxByBzB to oExEyEzE and oExEyEzE to oIxIyIzI can be expressed, respectively, in oBxByBzB as 2 3 2 3 2 3 0 0 0 6 7 6 7 6 7 T T T B xBE ¼ Rz ðnÞRy ðaÞ4 0 5 þ Rz ðnÞ4 a_ 5 þ 4 0 5 d_ 0 n_ 2 3 ð22Þ _ _ dsacn asn 6 7 _ ¼ 4 acn 5 _ þ dsasn _dca þ n_ and
2 3 2 39 0 > h_ = 6 7 6 7 6 7 T T T ¼ REB Ry ðp=2 þ uÞRx ðp=2 þ hÞ4 0 5 þ Ry ðp=2 þ uÞ4 0 5 þ 4 u_ 5 > > ; : 0 w_ 0 2 3 _ wcushÞ _ _ _ þ wsushÞðcdcn _ _ ðsdsn cacdcnÞðhsu ðcnsd þ cacdsnÞðu_ þ wchÞ þ cdsaðhcu casdsnÞðu_ þ wchÞ 6 7 _ wcushÞ _ _ þ wsushÞ _ ¼4 5 ðcdsn þ cacnsdÞðhsu þ sdsaðhcu 8 > <
2
0
3
_ þ wsushÞ _ _ _ wcushÞ _ caðhcu þ sasnðu_ þ wchÞ þ cnsaðhsu
ð23Þ
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
Matrices Rx, Ry and Rz are the rotating matrices along the x-, y- and z-axes, respectively. The superscript to the left of the symbol Bx means that the vector x is expressed in oBxByBzB. Denote the angular momentum of the sail as H = JBxBI, where the inertial moment matrix J = diag(Jx, Jy, Jz) is expressed in oBxByBzB and Jx, Jy and Jz and are the inertia moments about the xB, yB and zB axes, respectively. The attitude dynamics can be expressed using the following equation of the theorem of angular momentum B
_ þ B xBI H ¼ B T R þ eT TE ¼ H ¼ B CT qs þ B bT þ eT
s:t:qsi 2 f0; 1g
ð24Þ
where BTE represents the expected external torque expressed in oBxByBzB, and eT denotes the error between BTE and the torque generated from the RCDs. By substituting Eqs. (22) and (23) into Eq. (24), the equation can be rewritten in the form of Euler angles d, a and n. The unit vector along the yE axis and oINI expressed in oIxIyIzI is denoted as IyE and InI, respectively, which can be obtained as I
yE ¼ I r I v=kI r I vk2
ð25Þ
and I
nI ¼ ½ 0
0
ð26Þ
T
1 I yE
where vector Iv is the velocity of the sail expressed in oIxIyIzI. Furthermore, the orbital angles w, h and u can be obtained as 1 Þ
ð27Þ
0 0 T I nI Þ
ð28Þ
h ¼ arccosðI yE ½ 0 w ¼ arccosð½ 1
0
T
and u ¼ arccosðI nR I nI Þ
ð29Þ
By substituting BnR = REBRIEInR and Bn = REBRIEIn into Eq. (16), the solar torque BTR can be obtained and the theorem of angular momentum can finally be expressed in oBxByBzB. 3. Optimization for switch regulation of RCDs 3.1. Underdetermined Linear System with Binary Set Constraint The switch programming aims to attain an appropriate switch sequence for RCDs to achieve the desired control
qs
"
force and torque. Denote the programming vector Y = [msas, BTE]T 2 R61, and the equation of coupled attitude–orbit dynamics can be written as the following linear form Y ¼ Cqs þ b þ e
s:t:qsi 2 f0; 1g
ð30Þ
where matrix C, vector b and vector e can be represented as "I # "I # CF bF 6k C¼ 2R ; b¼ 2 R61 ; B B CT bT ð31Þ " # eF e¼ 2 R61 eT Naturally, the optimal value q*s of the power vector is sought to minimize the orbit and attitude errors, which can be expressed as the following multi-objective LSE (MLSE) problem 1 2 2 qs ¼ argqsi 2f0;1g min ðkeF k2 þ kT keT k2 Þ 2
ð32Þ
in which the denotation || ||2 represents the 2-norm of the vector, and kT is the Lagrange multiplier for weighting the attitude error. However, taking the electric energy consumption into account, it is expected to utilize as few switched-on RMDs as possible, if only the control accuracy is guaranteed; as a result, an l1–l2 hybrid optimization is developed as 1 kT 2 2 keF k2 þ keT k2 þ k1 kqs k1 qs ¼ argqsi 2f0;1g min ð33Þ 2 2 where the denotation || ||1 represents the 1-norm of the vector, and k1 is the Lagrange multiplier weighting the electric energy consumption. Furthermore, to provide the MLSE and the l1–l2 hybrid optimization with a unified quadratic programming (QP) scheme shown as follows, 1 T z Hz þ cT z z ¼ arg min zi 2f0;1g;16i6k 2 zj 2f1;0g;kþ16j62k
¼ arg
min
zi 2f0;1g;16i6k zj 2f1;0g;kþ16j62k
GðzÞ
ð34Þ
new matrices and vectors are defined as
ðCF þ kT CT ÞT ðCF þ kT CT Þ ðCF þ kT CT ÞT ðCF þ kT CT Þ
#
; ;H¼ T T qs ðCF þ kT CT Þ ðCF þ kT CT ÞðCF þ kT CT Þ ðCF þ kT CT Þ # " 1k1 CTF PF kT CTT PT c ¼ k1 þ 1k1 CTF PF þ kT CTT PT z¼
1151
ð35Þ
1152
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
Note that every element of power vector qs takes a value from the binary set {0, 1}, which is called the BSC. In addition, in most situations (when k P 6) the matrix C is underdetermined and does not have a Moore–Penrose inverse (a kind of pseudo inverse matrix, see (Zhang, 2004) for details). Therefore, the switch law design of RCDs is actually an underdetermined QP with BSC and has no analytical solution. As a result, a numerical method for the optimal search is sought. The problem has been converted into a QP in the previous discussion. Nevertheless, the BSC leads to an obstacle for directly applying the widely available algorithms (such as interior-point method, trust-region-reflectivity method, active-set method and so on) which rely on the first-order essential condition (i.e. 5G(z*) = 0) of the optimal solution, so the relaxation of the BSC to a continuous interval is performed before the optimal search. Meanwhile, a complete relaxation of BSC to the real number field will cause numerical solutions converging to an improper one. Thus, the BSC should be relaxed to the closed interval [1, +1] and an optimization problem can be derived from Eq. (34) as zs ¼ arg
min
zi 2½0;1;16i6k zj 2½1;0;kþ16j62k
GðzÞ
ð36Þ
in which the vector zs is a sub-optimal solution of problem (34). Furthermore, an approximate solution of problem (34) can be obtained after the discretization of zsi to the binary set {0, 1}, which can be written as ( 0 zsi < 0:5 qs ¼ ½qsi k1 ; qsi ¼ ; i ¼ 1; . . . ; k : 1 zsi P 0:5 ð37Þ
Fig. 3. The procedure of the BGP method.
3.2. Basic gradient projection: the BGP method A type of basic gradient projection (BGP) method (Bertsekas, 1999; Ma´rio et al., 2007) is proposed for the numerical calculation of Eq. (36). The geometry illustration of the BGP in a 2-dimensional case is given by Fig. 2. The BGP achieves the optimal search along the gradient descend of function G(z) with the complete relaxation of the BSC to the region of real number. After every step of descend, the sub-optimal solution is projected to the finite interval [1, +1] and the continuous optimal solution zs can be obtained. Subsequently, make the discretization of zsi to the binary set {0, 1} and finally obtain the optimal solution q*s . Denote the projection operator Proj as Proj½v ¼ ½vf ¼ ½vi k1 ;
Fig. 2. The geometry illustration of the BGP in a 2-dimensional case.
vi ¼ midð0; vi ; 1Þ
ð38Þ
in which function mid(a, b, c) denotes the middle value of its three scalar arguments. Define r as the damping factor of the backtracking, h as the backtracking step of argument, and l as the backtracking step of function value. The procedure of the BGP method is depicted in Fig. 3 and can be presented as follows. k1 Step 0 (initialization): Choose q(0) | s 2 W = {qs 2 R qsi 2 [1, +1]}; r 2 (0, 1), l 2 (0, 1/2) and [hmin, hmax] satisfying 0 < hmin. For p = 0, 1, 2, . . . Step 1: Compute h(0): h(0) = ((g(p))Tg(p))/((Hg(p))THg(p)). Replace h(0) by mid(hmin, h(0), hmax).
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
1153
Table 1 Parameters of the sail. Parameter
Specification
Mass (ms) Sail radius (rs) Payload region radius (rl) Interior SRP radius (rpin) Exterior SRP radius (rpex) Interior unit surface area (DSin) Exterior unit surface area (DSex) Moments of inertia (Jx, Jy, Jz) Spin rate (X)
64 kg 30 m 1m 9.3 m 29.5 m 77 m2 1.5 m2 (14,400, 14,400, 28,800) kg m2 2.5 rpm
Table 2 Parameters of the Venus rendezvous mission.
Fig. 4. The structure of solar sail demo.
Step 2 (backtracking line search): Choose h(p) = h(0), (0) rh , r2h(0), . . . such that (p) (p) (p) (p) T (p) G([q(p) s h 5G(qs )]f) 6 G(qs ) l[5G(qs )] {qs (p) (p+1) (p) (p) [q(p) = [q(p) 5G s h 5G(qs )]f} and set qs s h (q(p) s )]f. Step 3 (termination test): End the search if the termination condition is satisfied; otherwise, return to Step 1 for a new iteration. Here, the termination condition is chosen as kGðqsðpþ1Þ Þ GðqðpÞ s Þk2 6 tolG
ð39Þ
Vector g is defined as ( ðpÞ ðrGðqðpÞ qsi > 0 or ðrGðqðpÞ ðpÞ s ÞÞi ; s ÞÞi < 0 gi ¼ 0; otherwise: (p)
ð40Þ 4. Numerical demonstration The coupled attitude–orbit dynamics of a solar sail utilizing RCDs is numerically calculated in a Venus rendezvous mission. This section will present and analyse numerical results of both MLSE and l1–l2 hybrid optimization given by Eqs. (32) and (33). 4.1. Sail and mission description A circular solar sail depicted as Fig. 4 is used as a demo for detailed numerical calculations. The sail is a rigid even circular sheet with radius rs and total mass ms (payload mass included). On the surface of the sail, there are 16 independently switched RCDs in which 8 units are fixed on the interior region where rpin is the radius from oB to the interior SRP center and DSin is the surface area for every interior device; meanwhile, the other 8 units are fixed on the
Parameter
Specification
Initial time (TI) Initial position (xI, yI, zI) Initial velocity (vxI, vyI, vzI) Terminal time (TT) Terminal position (xT, yT, zT) Terminal velocity (vxT, vyT, vzT) Transfer time (tT)
24 May 2016 04:19:03.360 UTCG (0.4622, 0.9032, 1.4647 106) AU (2.6008 104, 1.3591 104, 2.2335) m/s 12 Sep 2017 05:07:17.76 UTCG (0.0874, 0.7196, 0.0149) AU (3.4599 104, 4.4822 103, 1.9295 103) m/s 476.0335 day
edge of the sail where rpex is the radius from oB to the exterior SRP center and DSex is the surface area for every interior device. Surrounded by the interior RCD, a circular region with radius rl is designed to have payloads (e.g., satellite computer, camera, etc.) fastened on. The remaining region of the sail is mainly membrane to hold the structure and collect solar energy. Table 1 presents detailed values for parameters of the sail. Notice that DSi is much larger than DSe, while rpi is much smaller than rpe. Thus, the solar acceleration is mainly generated by interior RCDs, while the torque is mainly generated by the exterior RCDs. The attitude control and the orbit control can be partly decoupled in structure. Besides, the circular sheet configuration is beneficial for the on-orbital deployment and the location of payloads will result in small centrifugal force or turbulent torque. For the Venus rendezvous, it is assumed that the solar sail departs from the Earth with no escape velocity, i.e., C3 = 0 and no propulsion occurs at the end of the mission. The reference orbit and attitude are calculated using an indirect method (Gong and Li, 2014; Jayaraman, 1980; Wood, 1982) to optimize the transfer time, assuming that the reflectivity modulation ratio is continuous. Parameters of the Venus rendezvous mission are presented in Table 2. 4.2. Numerical results and analysis Consider the MLSE problem given by Eq. (33) to ensure the accuracy of orbit control as well as compensate some
1154
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
Terminal Position Earth zI ˄AU˅
Reference Orbit Practical Orbit
Initial Position
Venus
xI˄AU˅ yI˄AU˅
Fig. 5. Practical transfer orbit of the sail when kT = 0.5.
5
2
x 10
100
1
0
50
0
50
100
150
200
250 300 Time(day)
350
400
450
Velocity Error(m/s)
Position Error(km)
Position Error Velocity Error
0 500
Fig. 6. Position and velocity error of the transfer orbit when kT = 0.5.
part of desired torque. Set the Lagrange multiplier kT = 0.5 for weighting attitude error. The numerical results are illustrated as Figs. 5–10. Fig. 5 illustrates the reference and practical transfer orbits which coincide with each other in the figure and can hardly be distinguished. Therefore, it is necessary to evaluate the ability of the switch programming by Fig. 6 which illustrates increasing errors in position and velocity that finally reach approximately 105 km and 50 m/s, respectively. The order of the terminal position error indicates that the sail finally locates beneath the influence surface of the Venus; meanwhile the terminal velocity of the sail is close to the Venus. Thus, the numerical result is compatible with the reference mission. These errors can be reduced and finally eliminated by increasing the total amount and distributed surface area of RCDs and decreasing the surface area of every single RCD, as indicated by Fig. 14 and Mu et al. (2013, 2015a,b). Fig. 7 shows the expected solar acceleration and practical one generating from RCDs along three axes of the inertial coordinate. The SRP
Fig. 7. SRP acceleration and error when kT = 0.5.
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
1155
Fig. 8. SRP torque compensation when kT = 0.5.
acceleration is approximately 104 m/s2, while the error between the expected and practical SRP acceleration is approximately 107 m/s2. Both of the expected and practical SRP torques along three axes of the body-fixed coordinate are illustrated in Fig. 8, which indicates a partial compensation of the torque by utilizing RCDs. However, as the SRP force given by Eq. (4) has a component along vector nR, the component of the SRP torque along zB axis is approximately 103 Nm, which performs as a disturbance of the spinning stability. Figs. 9 and 10 illustrate the discrete switch law for the eight interior RCDs and the eight exterior RCDs, respectively. Clearly, the interior RCDs mainly contribute to the orbit control and are switched on for nearly the entire voyage; while the exterior RCDs mainly contribute to the attitude control and are switched between the on and off states much more
frequently. Furthermore, as the Lagrange multiplier kT varies, the optimization objective and numerical consequences of switch programming change, which are illustrated in Figs. 11–13. Define the maximum relative error of SRP acceleration as ea max ¼ maxðkeF k2 =kms as k2 Þ 100%
ð41Þ
where function max(s) represents the maximum value of the scalar s. Illustrated by Figs. 11 and 12, as kT increases from 0 to 1, eamax increases from approximately 0.03% to 0.055%; while terminal position and velocity error increases from approximately 1.5 105 km to 1.65 105 km and from approximately 48 m/s to 53 m/s, respectively. On the other hand, define the uncompensated time of SRP torque as
Fig. 9. Discrete switch law for interior RCDs when kT = 0.5.
1156
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158 1
Unit 1 0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0 1
50
100
150
200
250
300
350
400
450
500
0 0
50
100
150
200
250
300
350
400
450
500
Unit 2 Unit 3 Unit 4 Unit 5 Unit 6 Unit 7 Unit 8
Time (day)
Fig. 10. Discrete switch law for exterior RCDs when kT = 0.5.
0.055
0.05
400 350
0.045
Uncompensated time(day)
Maximum relative error of acceleration(%)
450
0.04
0.035
0.03
300 250 200 150 100
0.025
0
0.1
0.2
0.3
0.4
0.5 0.6 Lambda
0.7
0.8
0.9
1
50
0
0.1
0.2
0.3
0.4
Fig. 11. Maximum relative error of SRP acceleration versus kT.
0.5 0.6 Lambda
0.7
0.8
0.9
1
Fig. 13. Uncompensated time of SRP torque utilizing RCD versus kT.
5.5
x 10
Terminal Position Error Terminal Velocity Error
TerminalPosition Error Terminal Velocity Error
1.6
52
1.5
50
Position Error(log(km))
5
Velocity Error(m/s)
Position Error(km)
50
54
45 40 35
4.5 30 4
25 20
3.5
Velocity Error(m/s)
5
1.7
15 10
3 5
1.4 0
0.1
0.2
0.3
0.4
0.5 0.6 Lambda
0.7
0.8
0.9
48 1
Fig. 12. Terminal position and velocity error of rendezvous versus kT.
2.5 10
20
30
40
50
60
70
0 80
k
Fig. 14. Terminal position and velocity error of the Venus rendezvous versus k.
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
1157 80
1500
80
L1=0,Position Error L1=0,Switched-on RMD
1000
Postion Error(km) 500
79 Switched-on 78 RMD
79
78
00 00
50 50
100 100
150 150
250 250
200 200 200
300 300 300
350 350 350
400 400 400
450 450
10000
77 77 500 500 500 80
80
Postion Error(km) 5000
Switched-on 78 RMD
78
L1=5e-6,Position Error L1=5e-6,Switched-on RMD 00 00
50 50 50
100 100 100
150 150 150
250 250 250
200
200 200
300 300 300
350
350 350
400 400 400
450
450 450
76 76 500 500 500
5
2
Postion Error(km)
x 10
100
100
L1=1e-5,Position Error L1=1e-5,Switched-on RMD
Switched-on 50 RMD
1
00 00
50
50 50
100 100
150 150
200 200 200
250 250 Time (day)
300 300 300
350 350 350
400 400 400
450 450
0 0 500 500 500
Fig. 15. Position error and number of switched-on RCD unit for l1–l2 hybrid optimization.
tunc ¼ timeðkeT k2 =kB T E k2 > eT Þ
ð42Þ
where function time(S) represents the total time that the statement S is true, and eT represents the threshold of the truth of statement S. Here, eT is set to equal 0.5. It is illustrated in Fig. 13 that the uncompensated time of SRP torque decreases from approximately 400 days to 100 days when kT increases from 0 to 1. Numerical results in Figs. 11–13 clearly demonstrate the weighting effect of the Lagrange multiplier for MLSE. As the value of kT increases, the optimal switch law of RCDs strengthens the compensation of desired torque and weakens the accuracy of orbit control, which can also be deduced from the mathematical expression of MLSE shown as Eq. (32). Since the reference orbit and attitude are obtained by assuming that the modulation surface area of RCDs is continuous, a decrease of both terminal position and velocity errors is expected with an increase in the total amount of RCDs. This inference can be verified numerically in the Venus rendezvous mission. Assuming kT = 0 which means the MLSE only considers the orbital control accuracy, recycle the switch programming of RCDs for different k values (the total amount of RCDs) and calculate the terminal position and velocity errors. The numerical result is shown in Fig. 14 in which the left vertical axis ‘Position Error (log(km))’ represents the log10 value of the terminal position error. It is clear that versus the increase of k from 16 to 80, the terminal position error falls from approximately 10 km to 103 km and the terminal velocity error falls from approximately 50 m/s to no more than 1 m/s. Thus, as stated previously, it is reasonable to infer that the orbital control error will finally be eliminated by RCDs as the modulation ratio tends towards a continuous variable. Finally, the l1–l2 hybrid optimization defined as Eq. (33) is computed numerically. Set k = 80 and kT = 0.01 and consider three situations: k1 = 0, k1 = 5 106 and k1 = 105. The position error and the number of switched-on
RCDs are illustrated in Fig. 15. The increase of the Lagrange multiplier k1 adds the weight of the consumption of electrical energy for the optimization. This effect is clearly illustrated in the numerical results as an obvious increase of terminal position error from approximately 1500 km to 105 km as well as a decrease of amount of switched-on RCDs. 5. Conclusions This paper derived coupled attitude–orbit dynamics of a solar sail utilizing independently switched reflectivity control devices (RCDs). According to the derivation, the solar force and torque were linearly dependent on the switch vector qs which usually had a higher dimension than the summing dimension of the force and torque. In addition, every element of vector qs took a value from a binary set {0, 1} due to the switch property of RCDs. Thus, the coupled attitude–orbit dynamics could be described by an underdetermined linear system (ULS) subject to the binary set constraint (BSC) in mathematics. This discovery helped the authors formulate the switch programming of RCDs into a multi-objective least-squares estimation (MLSE) and an l1–l2 hybrid optimization, respectively. By reformulating these two optimization problems into a unified quadratic programming (QP), a basic gradient projection (BGP) method was designed for the search of optimal solutions. A circular rigid sail with RMDs flying in the Venus rendezvous mission was taken as a demonstration of the derived dynamics and proposed BGP. Numerical results indicated that with an appropriate switch law of RMDs, the sail achieved to arrive beneath the influence surface of the Venus with approximately 105 km position error and 50 m/s velocity error. RMDs also generated solar torque which could compensate some part of desired angular acceleration of the sail. Furthermore, numerical results also illustrated that the terminal position and velocity error decreased below approximately 103 km and 0.5 m/s,
1158
T. Hu et al. / Advances in Space Research 57 (2016) 1147–1158
respectively, as the total amount of RCDs increased to 80. This phenomenon indicated a promising way to reduce errors for orbital control of the sail. The switch programming of RCDs in this paper was a kind of open-loop control, so a feedback controller was expected to eliminate the errors in the future. On the other hand, every element of the switch vector took a value from the binary set {0, 1}, which was a simplification for theoretical analysis. As a matter of fact, values of optical coefficients varied as the wavelength, and accurate synthetic values was looked forward to by on-orbital experiments. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grants no. 11272004). The authors wish to thank Fangming Han from Tsinghua University for his introduction to the underdetermined linear system with a binary set constraint. They also appreciate Pascal Willis and both anonymous reviewers for their comments and kind help. References Aliasi, G., Mengali, G., Quarta, A., 2013. Artificial Lagrange points for solar sail with electrochromic material panels. J. Guidance Control Dyn. 36, 1544–1550. Bertsekas, D., 1999. Nonlinear Programming, Athena. Gong, S., Baoyin, H., Li, J., 2009. Coupled attitude–orbit dynamics and control for displaced solar orbits. Acta Astronaut. 65, 730–737. Gong, S., Li, J., 2014. Fuel consumption for interplanetary missions of solar sailing. Sci. China-Phys. Mech. Astron. 57, 521–531. Jayaraman, T., 1980. Time-optimal orbit transfer trajectory for solar sail spacecraft. J. Guidance Control Dyn. 3, 536–542. Macdonald, M., McInnes, C., 2011. Solar sail science mission applications and advancement. Adv. Space Res. 48, 1702–1716. Ma´rio, A., Robert, D., Stephen, J., 2007. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems, IEEE. J. Sel. Top. Signal Process. 1, 586–597.
McInnes, C.R., 1999. Solar Sailing: Technology, Dynamics and Mission Applications. Springer-Verlag, Berlin (pp. 90–95). Mu, J., Gong, S., Li, J., 2013. Reflectivity-controlled solar sail formation flying for magneto sphere mission. Aerosp. Sci. Technol. 30, 339–348. Mu, J., Gong, S., Li, J., 2015a. Coupled control of reflectivity modulated solar sail for GeoSail formation flying. J. Guidance Control Dyn. 38, 740–751. Mu, J., Gong, S., Ma, P., Li, J., 2015b. Dynamics and control of flexible spinning solar sails under reflectivity modulation. Adv. Space Res. 56, 1731–1751. Rios-Reyes, L., Scheeres, D., 2005. Generalized model for solar sails. J. Spacecraft Rockets 42, 182–185. Roxburgh, I., 2001. Gravitational multipole moments of the Sun determined from helioseismic estimates of the internal structure and rotation. A. & A. 377, 688–690. Scholz, C., Romagnoli, D., Dachwald, B., Theil, S., 2011. Performance analysis of an attitude control system for solar sails using sliding masses. Adv. Space Res. 48, 1822–1835. Tsuda, Y., Mori, O., Funase, R., Sawada, H., et al., 2011a. Flight status of IKAROS deep space solar sail demonstrator. Acta Astronaut. 69, 833– 840. Tsuda, Y., Saiki, T., Funase, R., Mimasu, Y., 2011b. On-orbit verification of fuel-free attitude control system for spinning solar sail utilizing solar radiation pressure. Adv. Space Res. 48, 1740–1746. Tsuda, Y., Saiki, T., Funase, R., Mimasu, Y., 2013. Generalized attitude model for spinning solar sail spacecraft. J. Guidance Control Dyn. 36, 967–974. Wie, B., 2004a. Solar sail attitude control and dynamics, Part 1. J. Guidance Control Dyn. 27, 526–535. Wie, B., 2004b. Solar sail attitude control and dynamics, Part 2. J. Guidance Control Dyn. 27, 536–544. Wood, L., Bauer, T., Zondervan, K., 1982. Comment on ‘‘time-optimal orbit transfer trajectory for solar sail spacecraft”. J. Guidance Control Dyn. 5, 221–224. Zeng, X., Alfriend, K.T., Vadali, S.R., 2014. Solar sail planar multireversal periodic orbits. J. Guidance Control Dyn. 37, 674–681. Zhang, J., Wang, T., 2013. Coupled attitude–orbit control of a flexible solar sail for a displaced solar orbit. J. Spacecraft Rockets 50, 675–685. Zhang, X., 2004. Matrix Analysis and Applications. Tsinghua University, Beijing (pp. 85–96).