New Sliding Mode Attitude Controller Design Based on Lumped Disturbance Bound Equation Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. A. Sofyalı, Ph.D. 1; and E. M. Jafarov, Ph.D., D.Sc. 2 Abstract: In this paper, first the mathematical models of four major environmental disturbance torque components and corresponding bound equations are presented. Then the effect of the inertia matrix uncertainty on rigid satellite’s attitude dynamics is defined as an external torque in the derived state equation, which together with environmental torques forms the so-called lumped disturbance torque. After obtaining the complete equation for the model uncertainty-induced torque vector, its exact bound equation is derived by using matrix-vector norm relations. To validate the significance of these preliminary results for use in robust attitude controller design, a new modification of the classical sliding mode attitude controller present in literature is proposed, which is the primary contribution of this paper. The new design that is based on comprehensive knowledge of the lumped disturbance’s bounded variation leads to a decision rule on the switching control gain that is not excessively conservative. After verifying the accuracy of the bound equations in a simulation under no control, a second simulation is carried out with control input from the designed sliding mode controller to show that the proposed design works. The superiority of the new design is discussed in comparison with another design from literature that does not exploit the complete model of the inertia matrix uncertainty-induced torque through a comparative simulation’s result. The conclusion is that the modified controller design results in an attitude control system that has guaranteed robust stability in addition to reasonable conservativeness thanks to the newly obtained comprehensive decision rule on the switching control gain. DOI: 10.1061/(ASCE)AS.1943-5525.0000791. © 2017 American Society of Civil Engineers. Author keywords: Environmental disturbance; Inertia matrix uncertainty; Lumped disturbance; Matrix-vector norm; Sliding mode attitude control. Introduction There are two main classes of uncertainty sources in satellite attitude dynamics: (1) environmental effects and (2) uncertainty in the information on moments and products of inertia. Controlling attitude under uncertain conditions requires either mathematical models of disturbing effects on attitude dynamics or information on their bounds if they are bounded. Although some robust control techniques such as sliding mode control does not necessitate modeling of unwanted effects, the maximum values they can take during the whole operation, namely the bounds, must be somehow known at the controller design step because, e.g., in the sliding mode control, the so-called reaching condition is satisfied only if the switching control gain can be selected to be higher than the bound of the lumped disturbance for the whole control process. If such information is absent or not reliable, which may be the case in the practical application of the sliding mode control method, the preferred solution is mostly the easier one, namely to decide on a gain value that is expected to be sufficiently high. This easier solution results in excessive conservativeness that leads to unnecessarily high control effort and rise in the undesired phenomenon called chattering, which means 1 Control System Design Engineer, EDS Aerospace, ARI-1 Teknokent, Istanbul Technical Univ., Ayazaga Campus, Maslak, Sariyer, Istanbul 34469, Turkey (corresponding author). E-mail: sofyali@hotmail.com 2 Professor, Faculty of Aeronautics and Astronautics, Dept. of Aeronautical Engineering, Istanbul Technical Univ., Ayazaga Campus, Maslak, Sariyer, Istanbul 34469, Turkey. Note. This manuscript was submitted on February 8, 2017; approved on May 26, 2017; published online on October 25, 2017. Discussion period open until March 25, 2018; separate discussions must be submitted for individual papers. This paper is part of the Journal of Aerospace Engineering, © ASCE, ISSN 0893-1321. © ASCE high-frequency oscillation in a system’s dynamic response to discontinuous control effect (Cong et al. 2012). Another solution is to estimate the unknown bound relying on an adaptation scheme, which is known in literature as the adaptive sliding mode control approach and is so far highly studied on (Hu et al. 2007). The adaptive solution may lead to boundless increase in the estimated gain value due to the fact that the sliding surface vector becomes not exactly equal to zero even if the sliding mode is entered, which is named as the parameter drift problem and dealt with by using the dead-zone technique, the σ-modification, or the boundary layer scheme (Cong et al. 2012; Hu et al. 2007). In addition, there is another problem called the over-adaptation emerging from overestimation of the gain due to, e.g., an unrelated reaching phase adaptation, which triggered another line of studies in literature (Cong et al. 2012). As shown, besides being more complex than the conventional sliding mode control, the adaptive version has also issues to be dealt with. This paper deals with uncertain attitude dynamics of a rigid satellite in a circular low-Earth orbit (LEO) under environmental disturbances and inertia matrix uncertainty by relying on the fact that those uncertainty effects are bounded. With the aim of relieving the need to estimate the switching control gain, this paper presents the infinity norm of the disturbance torque vector due to model uncertainty, which is calculated rigorously from the derived torque equation and the definitions for inertia matrix uncertainty by employing matrix-vector norm relations. The torque vector equation is derived as an explicit function of only the time variable and the states, whereas it is presented in literature as a function of also the angular velocity’s rate to the best knowledge of the authors (Chen et al. 2014; Zhao et al. 2012). Because equations of motions are not expressed in state space in Chen et al. (2014) and Zhao et al. (2012) or many others, the derivation of that vector is not complete. This work also proposes bound equations on all of four major environmental 04017082-1 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. disturbance torque components, namely the gravity-gradient, the aerodynamic drag, the solar pressure, and the residual magnetic torque vectors. Studies such as Yang and Sun (2002) also present the maximum possible values of external torque components, even though not of all the four; the contribution of this study is numerical verification of the given bounds’ validity through high-fidelity simulation. In simulations, exact or realistic mathematical models of all four external disturbance components are used. All of those four models are also included in the text, which is preferred by few papers such as Meng et al. (2012). The first purpose of this paper is to derive the bound of the lumped disturbance torque that consists of the external disturbance torque components and the undesirable torque due to satellite inertia matrix uncertainty. The bound of the lumped disturbance is equal to the sum of the infinity norms of the four environmental disturbance torque vectors and the satellite model uncertaintyinduced torque vector. Equations that give the infinity norms of the gravity-gradient, aerodynamic drag, solar pressure, and residual magnetic torques are provided, which require knowledge of the mission orbit and satellite’s mass and dimensional properties in addition to access to some well-known books on attitude dynamics and control. Relying on the relation for the inverse of a sum of matrices, the torque due to model uncertainty is derived as a function of the states, the control input, and the resultant external disturbance torque. The inertia matrix uncertainty is modeled by two additive matrices: the diagonal one corresponding only to the parametric uncertainty in principal moments of inertia; and the full one representing solely the difference due to the angular deviation of the body axis system from the principal axis system. By using those two definitions, the infinity norm of the model uncertainty-induced torque is calculated; the result is a function of two respective parametric uncertainty bounds, minimum and maximum principal moments of inertia, two norms of the angular velocity and control vectors, and infinity norms of the four environmental torque vectors. All derivations are carried out for both inertial and Earth (nadir) pointing cases. In addition, three different types of actuators are considered for the definitions. Sliding mode control is a powerful tool for providing the controlled system with robustness to parameter uncertainties and external disturbances, especially if the undesired effects are bounded and their bounds are known (Utkin 1992; Edwards and Spurgeon 1998; Jafarov 2009; Shtessel et al. 2014). Thus a new sliding mode attitude controller is designed by modifying a classical one from literature using the previously obtained results, which is the second and main purpose of this work. The resulting controller is compared to another one that is designed based on incomplete knowledge of the model uncertainty-induced torque, which is mostly the case in literature to the best knowledge of the authors. Simulations are carried out by employing exact or realistic mathematical models for torques and the Sun position vector together with high fidelity atmosphere and geomagnetic field models. The paper is organized as follows. First, the state equation of attitude motion under environmental disturbances and inertia matrix uncertainty is derived, which is manipulated to obtain a structure convenient for robust control application. The terms in the state equation are defined for both inertial and Earth pointing cases. Second, the mathematical models of four environmental torques and the equations for bound values of the signals in their three channels are presented. The bound equations are written regarding the fact that those torques are of infinity norm-bounded nature. Then, the infinity norm of the satellite model uncertainty torque is obtained by using matrix algebra and norm relations. In the following section, a new sliding mode attitude controller design is carried out. Finally, the results of three exemplary simulations are © ASCE presented. The first one is without control, the second one belongs to the designed control system, and in the third one, the result of a design approach that is common in literature is simulated for comparison purposes with the new controller. Equations of Attitude Motion under Disturbances and Parametric Model Uncertainty The attitude motion of a rigid satellite orbiting the Earth is described by a set of dynamic and kinematic differential equations, which are derived step-by-step in, for example, Wie (1998). To avoid trigonometric expressions leading to singularity, quaternions are preferred to attitude (Euler) angles to represent the angular orientation of the satellite body frame with respect to the reference frame. Then the orthogonal direction cosine matrix C between those two frames can be written as function of the vectorial [q ¼ ð q1 q2 q3 ÞT ] and the scalar (q4 ) quaternion components: 2 1 − 2ðq22 þ q23 Þ 6 C ¼ 4 2ðq1 q2 − q3 q4 Þ 2ðq1 q3 þ q2 q4 Þ 2ðq1 q2 þ q3 q4 Þ 2ðq1 q3 − q2 q4 Þ 3 1 − 2ðq21 þ q23 Þ 7 2ðq2 q3 þ q4 q1 Þ 5 ð1Þ 2ðq2 q3 − q4 q1 Þ 1 − 2ðq21 þ q22 Þ Dynamic and Kinematic Equations The three dynamic equations can be written as the following vectorial differential equation: J ω̇ þ ω × Jω ¼ Td þ Tc ð2Þ In this paper, J ≜ J n þ ΔJ is the 3 × 3 uncertain inertia matrix defined as the sum of the nominal inertia matrix J n and the inertia uncertainty matrix ΔJ; ω is the 3 × 1 angular velocity vector of the body reference frame B (Fig. 1) with respect to the inertial reference frame N (i.e., ECI) Td ¼ Tgg þ Taero þ Tsolar þ Tmag ð3Þ is the 3 × 1 environmental disturbance torque; and Tc is the 3 × 1 control torque. The four components of Td will be covered extensively. The vectorial 1 q̇ ¼ ðq4 ω þ q × ωÞ 2 ð4aÞ 1 q̇4 ¼ − ðq · ωÞ 2 ð4bÞ and the scalar differential equations describe the kinematics. Eqs. (4a) and (4b) are valid if the reference frame is N; in other words, inertial pointing (IP) is the case of analysis/synthesis. In that case, the direction cosine matrix is defined as C ≡ CB=N ≜ ½ n1 n2 n3 , where the vectors ni ði ¼ 1, 2; 3Þ are unit vectors of N written in B. If the reference frame is selected to be the orbital reference frame A, kinematic equations have to be rewritten as 1 q̇ ¼ ðq4 ωB=A þ q × ωB=A Þ 2 1 n ¼ ðq4 ω þ q × ωÞ þ ðq4 a2 þ q × a2 Þ 2 2 04017082-2 J. Aerosp. Eng., 2018, 31(1): 04017082 ð5aÞ J. Aerosp. Eng. 1 1 n q̇4 ¼ − ðq · ωB=A Þ ¼ − ðq · ωÞ − ðq · a2 Þ 2 2 2 ð5bÞ Case 1: Inertial Pointing The IP reference state is xN ≜ ½ 01×3 according to the following relation with n being the angular speed of the satellite around the Earth, namely the mean motion Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. ω ≡ ωB=N ¼ ωB=A þ ωA=N ¼ ωB=A − na2 ¼ ½fn ðxÞ þ ΔfðxÞ þ ½bn þ ΔbuðxÞ þ ½dn ðx; tÞ þ Δdðx; tÞ ð7Þ The 7 × 1 state vector is defined as x ¼ ½ qT q4 ωT T 2 6 ΔfðxÞ ¼ 6 4 03×3 7 7 5 0 −J −1 n ðω × ΔJωÞ − J −1 unc ðω ð13Þ × JωÞ is the 7 × 1 uncertain system vector 6 dn ðx; tÞ ¼ 6 4 2 6 ¼6 4 03×1 0 J −1 n Td ðx; tÞ 3 7 7 5 3 03×1 0 J −1 n ½Tgg ðxÞ þ Taero ðx; tÞ þ Tsolar ðx; tÞ þ Tmag ðx; tÞ 7 7 5 ð14Þ is the 7 × 1 nominal disturbance vector, and 2 6 Δdðx; tÞ ¼ 6 4 3 03×1 0 J −1 unc ½Tgg ðxÞ þ Taero ðx; tÞ þ Tsolar ðx; tÞ þ Tmag ðx; tÞ ð15Þ is the 7 × 1 uncertain disturbance vector. 3 7 6 bn ¼ 4 01×3 5 3 03×1 ð8Þ The 7 × 3 control matrix consists of the nominal 2 ð12Þ is the 7 × 1 nominal system vector 2 ẋðtÞ ¼ fðxÞ þ buðxÞ þ dðx; tÞ ð11Þ 3 1 6 2 ðq4 ω þ q × ωÞ 7 6 7 6 7 1 fn ðxÞ ¼ 6 7 6 − ðq · ωÞ 7 2 4 5 −J −1 n ðω × J n ωÞ Attitude State Equation with Parametric Uncertainty Representation of attitude equations in state space is more convenient to analyze the dynamical system. The properties such as being control-affine, time-varying, underactuated, being in regular form, and having matched/unmatched disturbances can be determined by studying the state equation. These properties also govern the controller synthesis by pointing out control methods that are suitable for the dynamical system. The state equation of disturbed attitude motion of a rigid satellite with uncertain inertia orbiting the Earth can be written as follows: 01×3 T 2 ð6Þ Then the considered problem is Earth pointing (EP) or nadir pointing with the corresponding direction cosine matrix C ≡ CB=A ≜ ½ a1 a2 a3 . Note 1. In this paper, the nominal inertia matrix is considered to be equivalent to J when B coincides with the principal reference frame. Thus J n ≜ diagðJ 1 ; J 2 ; J 3 Þ, where J i ði ¼ 1, 2; 3Þ is satellite’s ith principal moment of inertia. In reality, there is hardly such a coincidence. In addition, J i have parametric uncertainty due to imperfect computer-aided design (CAD) calculations or erroneous measurements even if fuel consumption/sloshing and/or mass movements are not present. To take these two kinds of uncertainties into account, ΔJ has to be a full matrix rather than a diagonal one. Its definition will be given together with the variation intervals of its uncertain elements based on reasonable assumptions. 1 ð9Þ Case 2: Earth (Nadir) Pointing The EP reference state is J −1 n xA ≜ ½ 01×3 1 0 −n 0 T ð16Þ and the uncertain 2 03×3 3 7 6 Δb ¼ 4 01×3 5 The definitions are as follows: ð10Þ J −1 unc parts. The derivation of J −1 unc seen in Eq. (10) will be provided. Other terms in Eq. (7) differ for IP and EP problems. Therefore their definitions will be given separately for those two cases. © ASCE 3 1 n ðq ðq ω þ q × ωÞ þ a þ q × a Þ 2 7 62 4 2 4 2 7 6 7 6 1 n fn ðxÞ ¼ 6 7 − ðq · ωÞ − ðq · a2 Þ 7 6 2 2 5 4 −1 2 −1 −J n ðω × J n ωÞ þ 3n J n ða3 × J n a3 Þ 04017082-3 J. Aerosp. Eng., 2018, 31(1): 04017082 2 ð17Þ J. Aerosp. Eng. 7 7 5 2 6 6 6 ΔfðxÞ ¼ 6 6 4 0 −J −1 n ðω ::: þ × ΔJωÞ − 3n2 ½J −1 n ða3 2 6 dn ðx; tÞ ¼ 6 4 Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. 3 03×1 × ΔJa3 Þ þ × JωÞ : : : J −1 unc ða3 × Ja3 Þ 0 J −1 n ½Taero ðx; tÞ þ Tsolar ðx; tÞ þ Tmag ðx; tÞ ð27Þ whereas Tunc ðx; tÞ ¼ −ω × ΔJω þ 3n2 ða3 × ΔJa3 Þ 7 7 5 0 þ Tsolar ðx; tÞ þ Tmag ðx; tÞ 7 7 5 ð20Þ Remark 1. The definitions in Eqs. (17)–(20) are valid only for a triinertial rigid satellite in a circular orbit with principal moments of inertia satisfying J2 > J1 > J3 ð21Þ because only then the gravity-gradient torque acts as a restoring torque about xA given in Eq. (16) (Hughes 2004). Otherwise, for both IP and EP problems, Eqs. (13)–(15) hold instead of Eqs. (18)–(20), and the last three rows in Eq. (17) are same with the ones in Eq. (12). Improved State Space Representation The aim of this subsection is to represent all terms in Eq. (7) that are functions of ΔJ with a single disturbance term, namely a 7 × 1 disturbance vector due to satellite model uncertainty, as dunc ðx; tÞ ≜ ΔfðxÞ þ ΔbuðxÞ þ Δdðx; tÞ ð22Þ As a result, the uncertain state equation takes the following form: ẋðtÞ ¼ fn ðxÞ þ bn uðxÞ þ dðx; tÞ ≜ fn ðxÞ þ bn uðxÞ þ dn ðx; tÞ þ dunc ðx; tÞ ð23Þ The right-hand side of Eq. (23) has as many terms as the one of the initial state equation in the first line of Eq. (7). After substitutions and by using the following expression for J −1 unc −1 −1 J −1 unc ¼ −J n ΔJJ − ΔJJ −1 ½−ω × Jω þ Tc þ Td ðx; tÞ ð19Þ 3 03×1 J −1 unc ½Taero ðx; tÞ Tunc ðx; tÞ ¼ −ω × ΔJω − ΔJJ−1 ½−ω × Jω þ Tc þ Td ðx; tÞ 3 03×1 2 6 Δdðx; tÞ ¼ 6 4 J −1 unc ðω 7 7 7 7 ð18Þ 7 5 ð24Þ ð28Þ is valid for the case of EP. Eq. (3) substitutes for Td ðx; tÞ in both cases. Remark 2. With derivation of Eq. (27) or Eq. (28), it is achieved to represent the effect of inertia uncertainty on attitude motion by a single, newly defined disturbance term, which can be simply added to Td ðx; tÞ in the right-hand side of Eq. (2) allowing J’s in its lefthand side replaced by J n ’s. The disturbance torque due to satellite model uncertainty is dependent on the angular velocity, the environmental disturbance torque, and also the control input. In attitude motion, the environmental disturbance torque components are infinity norm-bounded as will be shown in the next section, therefore Tunc ðx; tÞ is also expected to have the same norm property. Remark 3. Considering the given expressions for the terms in Eq. (23) together with its structure, it can be concluded that attitude control systems are control-affine and in regular form regardless of the control actuation type. Even though the definitions presented earlier are valid for attitude dynamics actuated by reaction thrusters (gas jets), the structure of Eq. (23) is preserved for systems with magnetic actuators or reaction wheels. The uncertain state equation of attitude control system with reaction wheels has the extra terms of −J −1 n ðω × hÞ and ΔJJ −1 ðω × hÞ in the equations of fn ðxÞ’s last three rows and Tunc ðx; tÞ, respectively. In this paper, h is the 3 × 1 resultant angular momentum vector stored in the reaction wheels placed along the three body axes (Wie 1998). Because there is no difference in definition of the control matrix, it can be stated that attitude control systems with either thrusters or wheels have matched disturbances. This property together with being control-affine and in regular f orm is required by robust control methods such as conventional (classical) or integral sliding mode control. If attitude dynamics are controlled by magnetic actuators via interaction with local geomagnetic field, the control matrix becomes b ¼ bðx; tÞ ¼ ½ 03×3 03×1 ðJ −1 CB ÞT T , where CB ¼ ~ tÞ=kBðx; tÞk22 with Bðx; tÞ being the 3 × 1 local geoB~ T ðx; tÞBðx; ~ magnetic field vector and Bðx; tÞ its skew-symmetric matrix. Such a difference renders the attitude control system time-varying and underactuated. In addition, disturbances are unmatched in the attitude control system with magnetic actuators (Sofyalı et al. 2015). which is derived from Henderson and Searle (1981) −1 J −1 ¼ ðJ n þ ΔJÞ−1 ¼ J −1 n ðI − ΔJJ Þ −1 −1 ¼ ðI − J −1 ΔJÞJ −1 n ≜ J n þ J unc Environmental Disturbances ð25Þ and dunc ðx; tÞ can be obtained in the same form with dn ðx; tÞ 3 2 03×1 7 6 7 ð26Þ 0 dunc ðx; tÞ ¼ 6 5 4 −1 J n Tunc ðx; tÞ In this paper, Tunc ðx; tÞ is the 3 × 1 disturbance torque due to satellite model uncertainty. When the case is IP © ASCE There are four major components of the environmental disturbance torque vector as shown in Eq. (3): gravity-gradient torque Tgg , aerodynamic drag torque Taero , solar pressure torque Tsolar , and residual magnetic torque Tmag . Gravity-Gradient Torque The gravity-gradient torque emerges from spatially varying interaction of satellite’s mass distributed through its volume with Earth’s gravity field. It has the following simple formula (Wie 1998): 04017082-4 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. Tgg ðxÞ ¼ 3n2 ða3 × J n a3 Þ ð29Þ Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. Remark 4. Vector a3 is the third column of CB=A . In the EP problem, C ≡ CB=A and a3 is simply the third column of C, which is given in Eq. (1). Conversely, CB=A has to be calculated in the IP problem due to the equivalence of C ≡ CB=N . For circular orbits, CB=N ½ v̂ ðv̂ × x̂Þ −x̂ is the necessary transformation formula, where x̂ and v̂ are respectively the unit vectors of satellite’s position and velocity vectors written in N. The mean motion for a circular orbit can be calculated by Eq. (30) sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 3.986 × 1014 rad nðhÞ ¼ ð6.378 × 106 þ hÞ3 s Fig. 1. Reference frames of attitude motion ð30Þ as a function of the orbital altitude h, which is in meters. The highest absolute value that one of the three components of Tgg , which is a vector signal, can get during the process, namely Tgg ‘s infinity norm is (Fortescue et al. 2011; Larson and Wertz 1999) gravity center along b1 , b2 , and b3 , respectively. In the utilized simulation environment, ρ is calculated by the high fidelity atmosphere model NRLMSISE-00 (Wertz 1978). The velocity vector in Eq. (32) is written in B, which is the result of the following matrixvector multiplication: vB ¼ CB=N v. The infinity norm of the vector signal Taero can be accepted to be o n kTgg ½xðtÞk∞ ¼ sup maxðjT gg ji ðτ ÞjÞ τ i 3 ¼ n2 maxðjJ 3 − J 2 j; jJ 3 − J 1 j; jJ 1 − J 2 jÞ; 2 i ¼ 1,2; 3 kTaero ½xðtÞ; tk∞ ¼ CD Aρv2 The aerodynamic drag torque is exerted to satellite’s body by atmospheric particles that impact on its ram surfaces during its orbital motion. This torque acting on satellites with shape of rectangular prism can be modeled by (Wertz 1978; Gregory 2004) 3 X The distance of the geometric center to the gravity center generally does not exceed one-fifth of the shortest edge length of satellites with shape of rectangular prism and with no retractable appendages. The CD is taken as equal to 2.5 both in Eqs. (32) and (33), which corresponds to the possibly worst case. The A in Eq. (33) is calculated as 50% higher than the product of two highest edge lengths. Satellite’s orbital velocity in a circular orbit is a function of h, which is in meters, as sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 3.986 × 1014 m vðhÞ ¼ 6.378 × 106 þ h s Tai ; i¼1 Ta1 Ta2 Ta3 82 9 3 1 > > > > sgnðv ÞX − ΔX >6 > 1 < = 7 2 1 6 7 B ; ¼ − CD YZρjv1 j 6 × v 7 −ΔY > > 4 5 2 > > > > : ; −ΔZ 82 9 3 −ΔX > > > > <6 1 = 7 1 B 6 7 sgnðv ÞY − ΔY ¼ − CD ZXρjv2 j 4 ×v ; 2 5 > > 2 2 > > : ; −ΔZ 82 9 3 −ΔX > > > > <6 = 7 1 B −ΔY 7×v ¼ − CD XYρjv3 j 6 4 5 > > 2 > > : 1 sgnðv ÞZ − ΔZ ; 3 2 ð34Þ and ρ in Eq. (33) is calculated by the exponential atmosphere model (Vallado 1997) h−h0 H ρðhÞ ¼ ρ0 e− ð32Þ where CD = aerodynamic drag coefficient of the satellite; X, Y, and Z = edge lengths along b1 , b2 , b3 (Fig. 1), respectively; ρ = local atmospheric density or the local total mass density as defined by the employed atmosphere model; vi = component of satellite’s velocity vector along bi ; and ΔX, ΔY, and ΔZ are the components of the position vector of satellite’s geometric center with respect to its © ASCE ð33Þ ð31Þ Aerodynamic Drag Torque Taero ¼ minðX; Y; ZÞ 5 ð35Þ which is valid for the altitude interval of 0 to 1,000 km. In this paper, ρ0 , h0 , and H are read for the orbital altitude from Table 7.4 in Vallado (1997). Solar Pressure Torque The solar pressure torque is applied to satellite’s body as a result of solar radiation pressure on its illuminated surfaces. This torque acting on satellites with shape of rectangular prism can be calculated by (Wertz 1978; Montenbruck and Gill 2000) 04017082-5 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. Tsolar ¼ 3 X T⊙i ; i¼1 Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. T⊙1 T⊙2 T⊙3 0 2 31 9 82 3 1 > > > > > 2 3 > 1 > 6 sgnBr̂B · 6 0 7CX − ΔX 7 > > 1 > @ 4 5 A > > 6 7 s=c⊙ <6 2 = 7 B 6 7 B 6 7 ¼ −CRP YZP⊙ r̂s=c⊙ · 4 0 5 × 6 × r̂s=c⊙ ; 0 7 > > > >6 7 > > > 4 5 −ΔY 0 > > > > > ; : −ΔZ 9 82 3 −ΔX > > > > 2 31 0 > 2 3 > > >6 7 0 > 0 > > > 6 7 = < 61 7 6 7 C B B 6 7 B B 7 · 1 Y − ΔY sgn r̂ 4 5 A @ ; ¼ −CRP ZXP⊙ r̂s=c⊙ · 4 1 5 × 6 × r̂ s=c⊙ s=c⊙ 6 7 > > > > 62 7 > > > > 5 0 0 > >4 > > ; : −ΔZ 9 82 3 −ΔX > > > > > > 2 3 >6 > 7 > 0 > −ΔY > 7 > = <6 6 7 0 2 3 1 B 6 7 B 6 7 0 ¼ −CRP XYP⊙ r̂s=c⊙ · 4 0 5 × 6 × r̂ s=c⊙ > 7 > >6 1 B 6 7C > 7 > > > 4 sgn@r̂Bs=c⊙ · 4 0 5AZ − ΔZ 5 1 > > > > > 2 ; : 1 which is derived in a similar fashion to the derivation of Taero in Gregory (2004). In this paper, CRP is the radiation pressure coefficient of the satellite; P⊙ ¼ 4.56 × 10−6 N=m2 is the solar radiation pressure value at 1 astronomical unit (AU) distance from the Sun (Montenbruck and Gill 2000); r̂Bs=c⊙ ¼ CB=N r̂s=c⊙ is the unit position vector of the Sun with respect to the satellite obtained after the subtraction r̂s=c⊙ ¼ r̂⊙ − x̂ is carried out where r̂⊙ is the unit position vector of the Sun with respect to the Earth and is calculated by the algorithm available in Vallado (1997). In this case Tsolar ’s infinity norm can be obtained by (Montenbruck and Gill 2000) kTsolar ½x; tk∞ 1.496 2 minðX; Y; ZÞ ¼ CRP AP⊙ 1.470 5 ð37Þ kTmag ½xðtÞ; tk∞ ¼ M res ð36Þ 2 × 8.1 × 1015 Tm3 ð6.378 × 106 þ hÞ3 ð39Þ Both in Eqs. (38) and (39), M res ¼ m × 1 × 10−3 A m2 =kg is used, which is the relation read from Table 8.1 in Chobotov (1991) regarding the fact that residual magnetic torque may dominate the other three environmental torques in small satellite missions, which is the dominant mission type in LEO. In this paper, m is satellite’s mass. Remark 5. If the mission orbit and the satellite’s dimensions and mass together with its nominal inertia matrix are known a priori, the bounds on environmental disturbance torques can be evaluated by Eqs. (31), (33), (37), and (39) during the design process of a robust control system provided that the utilized robust control approach does not necessitate instantaneous knowledge of disturbances during the controlled process. The variable CRP is taken as equal to 2, which is its maximum possible value, both in Eqs. (36) and (37). Again, A in Eq. (37) is 150% of the product of two highest edge lengths. Satellite Model with Parametric Uncertainties Residual Magnetic Torque Modeling The residual magnetic torque emerges from the interaction of the residual magnetic dipole moment M res of the satellite with the local geomagnetic field B according to Definition. In this paper, the inertia uncertainty matrix ΔJ is defined as follows: ΔJ ≜ ΔJ 1 þ ΔJ 2 Tmag ¼ Mres × B ð38Þ where both vectors are written in the body reference frame. The M res is inducted in the satellite body mainly by current loops. The B is computed in the simulation environment by the high fidelity spherical harmonic geomagnetic field model IGRF2015 (Davis 2004). The infinity norm of Tmag is the result of the multiplication of the local geomagnetic field’s magnitude calculated according to the simple dipole model of the field (Sidi 1997) with M res © ASCE ð40Þ The first term 2 δ1 J1 6 ΔJ 1 ≜ 4 0 0 0 δ2 J2 0 0 3 7 0 5 ð41Þ δ3 J3 represents the parametric uncertainty due to imperfect CAD calculations or erroneous measurements of principal moments of inertia of a rigid satellite. The second term 04017082-6 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. 2 δ11 6 ΔJ 2 ≜ ð1 þ δ̄ J1 ÞmaxðJ 1 ; J 2 ; J 3 Þ4 δ12 δ13 δ 12 δ 13 3 δ 22 7 δ 23 5 δ 23 δ 33 ðkTgg k2 þ kTaero k2 þ kTsolar k2 þ kTmag k2 Þ pﬃﬃﬃ ≤ 3ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ þ kTmag k∞ Þ ð42Þ ð45Þ Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. From matrix algebra models the difference between the actual principal inertia matrix and the actual inertia matrix that is calculated with respect to the body reference frame B by six uncertain parameters. Assumption 1. The actual principal moments of inertia differ from their corresponding nominal values at most by 10%: δ 1 ; δ 2 ; δ 3 ∈ ½−0.1; þ0.1 ≜ ½−δ̄ J1 ; þδ̄ J1 . Assumption 2. The absolute value of the elements of ΔJ 2 does not exceed 15% of the value of ð1 þ δ̄J1 ÞmaxðJ 1 ; J 2 ; J 3 Þ: δ ij ∈ ½−0.15; þ0.15 ≜ ½−δ̄J2 ; þδ̄J2 ;i; j ¼ 1,2; 3. To obtain this bound, the principal frame is assumed to coincide with B after rotated diagonal axis, which has the unity vector of pﬃﬃﬃ around pﬃﬃﬃ its p ﬃﬃﬃ ½ 1= 3 1= 3 1= 3 T , for 10 degrees. A conservative approach is preferred to count for uneven inertia distribution, which is the case with elongated (e.g., gravity-gradiently stabilized) or flat (e.g., spin-stabilized) satellites. For satellites with evenly distributed inertia, i.e., with moments of inertia values close to each other, a much lower δ̄J2 value such as 0.06 is acceptable. k − Jki2 ¼ kJki2 ¼ ð1 þ δ̄J1 ÞJ max kJ −1 ki2 ¼ 1 1 1 ≜ ≜ ð1 − δ̄ J1 ÞminðJ 1 ; J 2 ; J 3 Þ ð1 − δ̄ J1 ÞJ min L1 k − ΔJki2 ¼ kΔJki2 ¼ ½δ̄ J1 þ δ̄ J2 ð1 þ δ̄ J1 ÞmaxðJ 1 ; J 2 ; J 3 Þ ≜ ðδ̄ J1 þ δ̄J2 þ δ̄ J1 δ̄ J2 ÞJ max ≜ L2 ( L2 ½L1 þ ð1 þ δ̄J1 ÞJ max kωk22 kTunc ½x; u; tk∞ ¼ L1 ! ~ i2 k − Jki2 kωk2 þ kuk2 : : : kωk : : : þ kTgg k2 þ kTaero k2 þ kTsolar k2 þ kTmag k2 ð43Þ Knowing that the induced two norm of a vector’s skewsymmetric matrix is equal to the two norm of that vector, ~ i2 ¼ kωk2 , Eq. (43) can be rewritten as i.e., kωk kTunc k2 ≤ L2 L1 ½L1 þ ð1 þ δ̄ J1 ÞJ max kωk22 kTgg k∞ þ kTaero k∞ : : : : : : þ kTsolar k∞ þ kTmag k∞ ! ) þ kuk2 þ kuk2 ð49Þ Case 2: Earth Pointing By starting with Eq. (28) ( L2 kTunc ½x; u; tk∞ ¼ ½L1 þ ð1 þ δ̄J1 ÞJ max kωk22 þ 3n2 L1 L1 ! ) kTgg k∞ þ kTaero k∞ : : : pﬃﬃﬃ þ kuk2 þ 3 : : : þ kTsolar k∞ þ kTmag k∞ is derived, which’s only difference from Eq. (49) is the presence of the extra term 3n2 L2 . Remark 6. Eq. (49) or Eq. (50) has to be computed onboard at each step of the control process to adjust the controller according to varying ω and u. The derivation results in this subsection are convenient for use in robust attitude controller design. Sliding Mode Attitude Controller Design 2 ð44Þ kq 6 sðxÞ ¼ ω þ 4 0 where the following norm inequality valid for vectors of size 3 × 1 is employed: © ASCE : : : þ kTsolar k∞ þ kTmag k∞ ) In this section, the design of a new sliding mode attitude controller will be carried out. The originality will be the usage of the derived infinity norms of the unwanted torque components in the design to obtain a decision rule on the switching control gain value that is not too conservative. The used sliding surface vector is ( pﬃﬃﬃ þ 3 ! kTgg k∞ þ kTaero k∞ : : : ð50Þ ~ i2 k − ΔJki2 kωk2 þ k − ΔJki2 kJ −1 ki2 kTunc k2 ≤ kωk × ð48Þ could be p derived or defined, which led to Eq. (44). The relation ﬃﬃﬃ kTunc k2 = 3 ≤ kTunc k∞ ≤ kTunc k2 holds between the two different norm types of 3 × 1 vectors, which allows taking the right-hand side of Eq. (44) as equal to kTunc k∞ pﬃﬃﬃ þ 3 Case 1: Inertial Pointing To derive the infinity norm of Tunc , an upper bound for its two (Euclid) norm is sought first because the induced two norms of the matrices in Eq. (27) are obtainable as follows: ð47Þ and Derivation of Bound on Disturbance Torque due to Inertia Matrix Uncertainty Note 2. In the following derivation, the norm inequality kTc k2 ≤ kuk2 is utilized. The equality is valid for control systems with gas jets or reaction wheels whereas the inequality holds for the torque Tc ¼ CB u induced by magnetic actuators because of the fact that the induced two norm of CB is equal to 1 (Sofyalı et al. 2015): kCB ki2 ¼ 1. ð46Þ 04017082-7 J. Aerosp. Eng., 2018, 31(1): 04017082 0 kq 0 3 07 5q; kq > 0 ð51Þ 0 0 kq |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} ≜K q J. Aerosp. Eng. which has become the preferred sliding vector in literature on attitude control after first proposed in Vadali (1986), where its structure is proved to guarantee that the sliding motion ends at the reference. In this case kq is the sliding surface design parameter. The utilized design approach is the equivalent control method, which results in a control vector that is the sum of the equivalent control vector ueq and the discontinuous reaching control vector ureach Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. Tc ≡ u ≜ ueq þ ureach ¼ sT ½ðTd þ Tunc Þ þ ureach ð52Þ The following reaching dynamics are assigned (Hung et al. 1993) as in Vadali (1986) and Lo and Chen (1995) to ureach ðxÞ ¼ −K ss ðxÞsgn½sðxÞ − K s sðxÞ; V̇ðsÞ ¼ sT J n ṡ ¼ sT J n ½Gðfn þ dn þ dunc Þ þ J −1 n u 2 3 1 J K ðq ω þ q × ωÞ − ω × J ω : : : n q 4 n 5 ¼ sT 4 2 : : : þ ðTd þ Tunc Þ þ ueq þ ureach 2 3 1 J K ðq ω þ q × ωÞ − ω × J ω þ ðT þ T Þ : : : n q 4 n d unc 5 ¼ sT 4 2 1 : : : − 2 J n K q ðq4 ω þ q × ωÞ þ ω × Jn ω þ ureach which becomes V̇ðsÞ ¼ sT ½ðTd þ Tunc Þ þ ðureach þ T unc Þ ¼ sT ½ðTd þ Tunc Þ þ ðI − ΔJJ −1 Þureach K ss ; K s > 0 ð53Þ where K ss ¼ diagðkss Þ and K s ¼ diagðks Þ are the continuous and discontinuous reaching control gain matrices of size 3 × 3, respectively, thus kss corresponds to the aforementioned switching control gain. Note also that the considered problem is IP, and the system is actuated by a reaction thruster triad, both the same with the cases in Vadali (1986). Derivation of the Equivalent Control Vector The equivalent control vector is derived under ideal sliding assumption, i.e., for s ¼ ṡ ¼ 03×1 and d ¼ 07×1 . If s is derived with respect to time by assuming that the attitude motion is in the sliding mode and no more subject to environmental disturbances and undesired effects due to model uncertainty ¼ −sT fðI − ΔJJ −1 Þ½K ss sgnðsÞ þ K s s − ðTd þ Tunc Þg ¼ −sT ½ðI − ΔJJ −1 ÞK ss sgnðsÞ − ðTd þ Tunc Þ − sT ðI − ΔJJ −1 ÞK s s Derivation of the Decision Rule on kss Eq. (59) indicates that kss > kðI − ΔJJ −1 Þ−1 ðTd þ Tunc Þk∞ is the first condition for the satisfaction of the reaching condition sT ṡ < 0. As known, for an arbitrary 3 × 1 vectoral signal a, the following inequalities hold between its infinity and two norms: pﬃﬃﬃ kak2 = 3 ≤ kak∞ ≤ kak2 ð60Þ Using Eq. (60) kðI − ΔJJ −1 Þ−1 ðTd þ Tunc Þk∞ ð54Þ ≤ kðI − ΔJJ −1 Þ−1 ðTd þ Tunc Þk2 can be obtained, where the 3 × 7 matrix G ≜ ∂s=∂x ¼ ½ K q 03×1 I and Gbn ¼ J −1 n according to Eq. (9). ueq can be solved for from Eq. (54) by using Eq. (12) as ≤ kðI − ΔJJ −1 Þ−1 ki2 ðkTd k2 þ kTunc k2 Þ L1 pﬃﬃﬃ 3ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ ≤ L1 − L2 þ kTmag k∞ Þ þ kTunc k2 ð55Þ In the reaching mode that precedes the sliding mode, s; ṡ ≠ 03×1 hold. If the ideality assumption is abandoned, the time derivative of s is derived as ṡðxÞ ¼ G½fn ðxÞ þ dðx; tÞ þ J −1 n uðxÞ ≠ 03×1 1 1 ≤ −1 1 − kΔJJ ki2 1 − kΔJki2 kJ −1 ki2 1 L1 ≤ ¼ 1 − L2 =L1 L1 − L2 kðI − ΔJJ −1 Þ−1 ki2 ≤ ð56Þ At this step, it is necessary to divide Tunc in Eq. (27) into continuous and discontinuous parts as follows: kTunc k2 ≤ ≜Tunc ð57Þ ≜T unc Differentiating the chosen Lyapunov function candidate VðsÞ ¼ ð1=2ÞsT J n s gives © ASCE ð62Þ Eq. (62) relies on the definitions in Eqs. (47) and (48), and implies that kΔJJ −1 ki2 ≤ L2 =L1 < 1 has always to be satisfied Tunc ¼ −ω × ΔJω − ΔJJ −1 ð−ω × Jω þ ueq þ Td Þ |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} þ ð−ΔJJ ureach Þ |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} ð61Þ can be obtained by using also the following inequality for the induced two norm of the inverse matrix ðI − ΔJJ −1 Þ−1 (Meyer 2000): Reaching Analysis −1 ð59Þ after the substitution of Eq. (57) into Eq. (58). ṡðxÞ ¼ Gẋ ¼ G½fn ðxÞ þ bn ðx; tÞueq ¼ Gfn ðxÞ þ J −1 n ueq ¼ 03×1 1 ueq ðxÞ ¼ − J n K q ðq4 ω þ q × ωÞ þ ω × J n ω 2 ð58Þ L2 f½L1 þ ð1 þ δ̄ J1 ÞJ max kωk22 L1 pﬃﬃﬃ þ 3ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ þ kTmag k∞ Þ þ kueq k2 g ð63Þ can be written in the same way with kTunc k2 by using its definition in Eq. (57). If the following inequality for ueq, which utilizes the ~ i2 ¼ 1 and kωk ~ i2 ¼ kωk2 equations kTðq; q4 Þki2 ≜ kq4 I þ qk 04017082-8 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. 1 ~ n ωk2 ≤ J max kωk22 kueq k2 ≤ k− J n K q Tðq; q4 Þωk þ kωJ 2 2 kq þ J max kωk2 ð64Þ 2 is used in Eq. (63), then the following inequality for kTunc k2 is arrived at: L kTunc k2 ≤ 2 ½L1 þ ð2 þ δ̄ J1 ÞJ max kωk22 L1 pﬃﬃﬃ þ 3ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ kq ð65Þ þ kTmag k∞ Þ þ J max kωk2 2 Finally, Eq. (65) is substituted into Eq. (61) to obtain kðI − ΔJJ −1 Þ−1 ðTd þ Tunc Þk∞ pﬃﬃﬃ L þ L2 ≤ 3 1 ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ þ kTmag k∞ Þ L1 − L2 kq L2 2 þ ½L1 þ ð2 þ δ̄ J1 ÞJ max kωk2 þ J max kωk2 L1 − L2 2 ð66Þ pﬃﬃﬃ L2 L2 kT unc k2 ≤ L1 ðkureach k2 Þ ≤ L1 ð 3kss þ ks ksk2 Þ, has to be added to the right-hand side of the inequality in Eq. (65). Then, that upper bound can be taken as equal to kTunc k∞ ¼ L2 ½L1 þ ð2 þ δ̄J1 ÞJ max kωk22 L1 pﬃﬃﬃ þ 3ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ þ kTmag k∞ Þ pﬃﬃﬃ kq þ J max kωk2 þ 3kss þ ks ksk2 ð69Þ 2 Derivation of the Decision Rule on ks The second reaching condition can be read from Eq. (59) as ðI − ΔJJ −1 ÞK s ≥ 0. It is straightforward to rewrite the condition as ðI − ΔJJ −1 ÞK s ¼ ½I − ðJ − J n ÞJ −1 K s ¼ ½ðI − IÞ þ J n J −1 K s ¼ J n J −1 K s ≥ 0 ð70Þ where J n , J −1 , and K s are symmetric and positive definite matrices by definition. Thus, J n J −1 K s is symmetric and positive definite. In conclusion, J n J −1 K s > 0 and the second reaching condition is satisfied for every value of ks . which renders the first reaching condition be written as kss ðtÞ ¼ kss ½xðtÞ pﬃﬃﬃ L þ L2 ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ þ kTmag k∞ Þ > 3 1 L1 − L2 kq L2 ½L1 þ ð2 þ δ̄J1 ÞJ max kωðtÞk22 þ J max kωðtÞk2 þ L1 − L2 2 ð67Þ Remark 7. The derivation of the rule in Eq. (67) is the primary contribution of this work, which is structured on complete modeling of the lumped disturbance and rigorous derivation of its bound. To the best knowledge of the authors, such a comprehensive mathematical relation between the switching control gain and the bounds on external disturbance components and inertia matrix uncertainty has not been presented yet in line of studies on sliding mode attitude control, which is started by the pioneering work of Vadali (1986) and followed by classical works such as Lo and Chen (1995). Remark 8. The rule in Eq. (67) is state dependent, so it enables adjusting kss according to the instant Euclid norm of the angular velocity vector during the operation as long it is available to the control system. This results in a reasonably conservative switching control action. Note 3. The following kss ðtÞ that satisfies the inequality in Eq. (67) can be used for the controller: pﬃﬃﬃ L þ L2 kss ðtÞ ¼ 1.1 × 3 1 ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ L1 − L2 L2 þ kTmag k∞ Þ þ ½L1 þ ð2 þ δ̄ J1 ÞJ max kωðtÞk22 L1 − L2 kq þ J max kωðtÞk2 ð68Þ 2 Note 4. To obtain the upper bound of the two norm of the model uncertainty-induced torque for the controlled case, the right-hand side of the inequality for kT unc k2 ‘s Euclid norm, which is equal to © ASCE Simulation Results and Superiority Discussion of the Proposed Design Simulation without Control The first exemplary simulation is carried out with u ¼ 0. The simulation orbit belongs to ITUpSAT-1, the picosatellite of Istanbul Technical University, which was launched on September 23, 2009. The model of the Danish microsatellite Ørsted is used. Table 1 presents the quantities belonging to the nearly circular simulation orbit (CelesTrak 2015) and Ørsted’s model (Wisniewski 1996). Note 5. The close values of its principal moments of inertia seen in Table 1 indicates that the used Ørsted model has an even inertia distribution, thus δ̄ J2 ¼ 0.06 is used in the evaluation of kTunc k∞ according to Assumption 2 together with δ̄ J1 ¼ 0.1 according to Assumption 1. The values of ΔX, ΔY, and ΔZ are chosen as 0.035, 0.025, and 0.05 m, respectively, which are lower than ð1=5Þ × 0.34 m. Table 1. Quantities Belonging to Simulation Orbit and Used Satellite Model Quantity Orbit h (km) n (degrees=s) v (km=s) ρ (kg=m3 ) Satellite model X (m) Y (m) Z (m) A (m2 ) m (kg) J 1 (kg m2 ) J 2 (kg m2 ) J 3 (kg m2 ) 04017082-9 J. Aerosp. Eng., 2018, 31(1): 04017082 Value 703.463 6.07 × 10−2 7.503 3.476 × 10−14 0.45 0.34 0.68 0.459 61.8 2.904 3.428 1.275 J. Aerosp. Eng. Table 2. Bounds on Disturbance Torques Distance torque components (Nm) 3,625 × 10−6 kTaero k∞ Eq. (33) 7,633 × 10−8 kTsolar k∞ Eq. (37) 2,948 × 10−7 kTmag k∞ Eq. (39) 2,819 × 10−6 kTunc ðtÞk∞ Eq. (49) Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. Value kTgg k∞ Eq. (31) 5,854 × 10−6 þ 2,439×kωðtÞk22 Fig. 3. Angular velocities Fig. 2. Quaternions Ørsted’s residual magnetic dipole moment vector is Mres ¼ ½ 0 3 × 10−2 0 T Am2 according to Bak (1999), which corresponds to an M res value lower than the one computed by M res ¼ m × 1 × 10−3 A m2 =kg. In the simulation, ΔJ 1 ≜ diagð0.1J 1 ; −0.1J 2 ; 0.05J 3 Þ is used, and ΔJ 2 is calculated by assuming that the angular difference between the body and the principal reference frames is 10 degrees around the diagonal axis of the principal frame. It is verified that Assumption 2 would not be violated even if Ørsted’s model with deployed gravity-gradient boom, which is ½J 1 ; J 2 ; J 3 ¼ ½181.25; 181.78; 1.28 kg m2 (Wisniewski 1996), would be used. The infinity norms of five disturbance torque components are tabulated in Table 2. The simulation is initiated at the reference state of the IP problem, which is given in Eq. (11). Figs. 2 and 3 show how the state variables are carried out of the reference state by the disturbance torques due to environmental effects and satellite model uncertainty during the first five orbital periods, each equal to 5,931 s or circa 99 min. The disturbance torques are shown in the following five figures as upper-bounded and lower-bounded by their infinity norm values. Remark 9. Only for demonstrational purposes, the gravitygradient torque in Fig. 4 is calculated by using J instead of J n : Tgg ðxÞ ¼ 3n2 ða3 × Ja3 Þ ð71Þ Thus the following related bound is superimposed on the graphs of Tgg ’s components. © ASCE Fig. 4. Gravity-gradient torque 3 kTgg k∞ ¼ n2 max½jð1 þ δ 3 ÞJ 3 − ð1 þ δ 2 ÞJ 2 j; jð1 þ δ3 ÞJ 3 2 − ð1 þ δ 1 ÞJ 1 j; jð1 þ δ 1 ÞJ 1 − ð1 þ δ2 ÞJ 2 j ð72Þ Figs. 2 and 3 point out that a control system that is robust against disturbances is highly necessitated to achieve IP or EP. Disturbances are observed to vary in the obtained bounds in Figs. 4–8, which means that the derivations in this paper are validated. The dependency of the varying bound in Fig. 8 on the two norms of the angular velocity is clear from a comparison between Figs. 3 and 8. The formulation of the uncertainty effect as a disturbance torque allows the comparison of its order with the orders of environmental disturbance torques. Under given simulation conditions, it is comparable to the dominant gravity-gradient torque and the residual magnetic torque. 04017082-10 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. Fig. 5. Aerodynamic drag torque Fig. 7. Residual magnetic torque Fig. 8. Model uncertainty torque Fig. 6. Solar pressure torque Table 3. Quantities Related to Simulation with Control and Initial Conditions Simulation with Control The second exemplary simulation is carried out with u ¼ ueq þ ureach . Table 3 presents the quantities related to the simulation. Note that kTgg k∞ , kTaero k∞ , kTsolar k∞ , and kTmag k∞ are the same in cases without and with control, which means that the mentioned four values in Table 1 are also valid for the second simulation. By using the same values for other quantities mentioned in the previous subsection, the following controlled attitude responses are obtained. As expected from a sliding mode controller, the disturbances due to environmental effects and model uncertainty are rejected (Fig. 9). In return, oscillations with high frequency and low amplitude are observable in angular velocity responses as shown in Fig. 10, which is induced by the discontinuous control signals. Fig. 11 depicts that the sliding mode is entered rapidly. Because the simulation step size is taken as 1 s, which corresponds to a © ASCE Quantity Control conditions δ̄J1 δ̄J2 L2 =L1 kq (1=s) ks (Nm s) kss ðtÞ (Nm) kTunc ðtÞk∞ (Nm) Initial conditions ðq1 ; q2 ; q3 ; q4 Þ0 ωT0 (×10−2 degrees=s) sT0 (×10−3 1=s) kω0 k2 (n) 04017082-11 J. Aerosp. Eng., 2018, 31(1): 04017082 Value 0.1 0.06 0.496 < 1 2.5 × 10−3 1 × 10−1 3.853 × 10−5 þ 8.211×kωðtÞk22 þ 4.215 × 10−3 ×kωðtÞk2 5.854 × 10−6 þ 8.211×kωðtÞk22 þ 4.215 × 10−3 × kωðtÞk2 þ 4.959 × 10−2 × ks½xðtÞk2 þ 8.589 × 10−1 × kss ðtÞ (0.123 0.707 0 0.696) (5 7 6) (1.18 2.99 1.047) 1.728 J. Aerosp. Eng. Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. Fig. 9. Controlled attitude angles Fig. 11. Sliding surface vector components Fig. 10. Controlled angular velocities Fig. 12. Model uncertainty torque under control switching frequency of 1 Hz that is a realizable value in application, the ideal sliding does not occur as observed in the third graph of Fig. 11. Ideal sliding requires infinite switching frequency, which is not realizable. In nonideal sliding the switching occurs not on the sliding manifold but in the some vicinity of the manifold. The higher the frequency, the narrower that vicinity. Fig. 12 shows that as angular velocity components converge to zero in steady state the bound on the model uncertainty-induced torque converges to their minimum values, which is constant as can be deduced from Table 3. Similarly, the switching control gain oscillates around a constant value in steady state (Fig. 13). Finally, the control torque components are presented in Fig. 14, which are in the order of 10−4 Nm. upper bounds on the disturbances due to environmental effects and inertia matrix uncertainty. In this subsection, a second rule on kss will be derived by using the Lyapunov function candidate VðsÞ ¼ ð1=2ÞsT Js instead of VðsÞ ¼ ð1=2ÞsT J n s together with the dynamic equations given in Eq. (2) as done in Chen and Lo (1993) and Lo and Chen (1995). Note that the term ΔJ ω̇ appears in the right-hand side of V̇ ’s equality if V ¼ ð1=2ÞsT J n s is employed, which renders the reaching analysis complicated. It is straightforward to obtain 2 6 C7 B 1 C7 − sT K s s B ~ V̇ ¼ −sT 6 4K ss sgnðsÞ − @Td þ 2 ΔJK q Tω − ωΔJω A5 |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} ≜TΔJ Superiority Evaluation through a Comparative Design and Simulation If the design would not be based on the state space model, the exact equation of Tunc , and its derived infinity norm, it might not be guaranteed that the derived decision rule on kss exceeds the total actual © ASCE 13 0 ð73Þ by using Eq. (2), Eq. (4a), the definition of the uncertain inertia matrix, Eqs. (53) and (55). Eq. (73) implies that the reaching 04017082-12 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. Fig. 13. Switching control gain Fig. 15. Continuous part of model uncertainty torque under control with Eq. (74) instead of Eq. (68) Fig. 14. Control torque condition is satisfied for kss > kTgg k∞ þ kTaero k∞ þ kTsolar k∞ þ kTmag k∞ þ kTΔJ k∞ . Following the necessary steps gives kss ðtÞ ¼ 1.1 × ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ þ kTmag k∞ Þ kq þ L2 kωðtÞk22 þ kωðtÞk2 2 ¼ 7.497 × 10−6 þ 5.691 × 10−1 ×kωðtÞk22 þ 7.113 × 10−4 × kωðtÞk2 ð74Þ which, regardless if ΔJ is zero matrix or not, reduces to kss ðtÞ ¼ 1.1 × ðkTgg k∞ þ kTaero k∞ þ kTsolar k∞ þ kTmag k∞ Þ in steady state, i.e., when ω → 0 as t → ∞. Thus the rule in Eq. (74) does not account for the disturbance due to inertia matrix uncertainty in steady state, which prevents the controller from producing discontinuous control action with adequate amplitude to counteract also the disturbance due to inertia matrix uncertainty besides the environmental ones. Conversely, the rule in Eq. (67) takes the nonzero © ASCE Fig. 16. Switching control gain with Eq. (74) instead of Eq. (68) Tunc ðx; ∞Þ in steady state into account thanks to the constant coefficient of ðL1 þ L2 Þ=ðL1 − L2 Þ, which is greater than one and becomes zero only if ΔJ is zero matrix. The same simulation as in the previous subsection is carried out by only using Eq. (74) instead of Eq. (68) for comparison purposes. The unwanted torque due to ΔJ is defined as TΔJ in Eq. (73). If even only the continuous part of Tunc in Eq. (27), namely Tunc according to Eq. (57), is plotted against the infinity norm of TΔJ , it can be numerically observed that kTΔJ ½xðtÞk∞ ¼ L2 ½kωðtÞk22 þ ðkq =2ÞkωðtÞk2 cannot serve as the bound equation on the model uncertainty-induced torque Tunc because kTΔJ ½xð∞Þk∞ ¼ 0 although ΔJ is not zero matrix in steady state (Fig. 15). Because the rule in Eq. (74) leads to a lower steady state value of kss compared to the one by the rule in Eq. (68) as clear from the comparison of Figs. 13–16, the steady state error margin of the attitude angles, which is approximately between −0.05° and þ0.05° 04017082-13 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. Conclusions Fig. 17. Controlled attitude angles with Eq. (74) instead of Eq. (68) In this work, a sliding mode attitude controller was designed on the basis of a newly derived decision rule on the discontinuous reaching law design parameter kss. In the derivation of the new rule, the lumped disturbance bound equation was utilized. To be able to combine the environmental disturbances and the model uncertainty effect into a lumped expression, the exact mathematical formula of the torque due to inertia matrix uncertainty was derived by using a novel complete uncertainty model. The primary contribution of this study is the structure of the aforementioned new decision rule, which mathematically guarantees that kss exceeds the maximum lumped disturbance torque component in any channel along the whole control process. This theoretically means that the reaching condition is satisfied and the sliding mode is preserved. It also renders kss variable according to the angular velocity vector and is thus not too conservative. Both inertial and Earth-pointing cases were covered throughout the text. Moreover, the differences in terms of the state equation in cases of employment of reaction thrusters, reaction wheels, or magnetic actuators were pointed out where necessary. A simulation under no control effect was carried out to verify the validity of all the derivations and assumptions presented in the text. Its results numerically verified that the bounds calculated by the related equations are neither exceeded by the corresponding vector signal components nor turned out to be excessively conservative. A second simulation numerically verified the functionality of the control system and the derived new rule. It showed that the classical sliding mode control method based on the equivalent control approach could be modified to robustly control the uncertain attitude dynamics in a reasonably conservative way thanks to the derived bound value of lumped disturbance, which is available to the robustly stabilizing controller at each instance when the angular velocity is known. Finally, the superiority of the new controller was clarified by comparing its decision rule on kss to the one of another controller that was designed without using the complete model of Tunc . A comparative simulation depicted the insufficiency of the other controller to produce discontinuous control signals with adequate amplitude in steady state. Acknowledgments Fig. 18. Control torque with Eq. (74) instead of Eq. (68) (Fig. 17), is slightly narrower than the margin obtained by the proposed controller and shown in Fig. 9, which is circa ½−0.2°; þ0.1°. Fig. 18 depicts that the initial magnitudes of the control torque’s components are lower by the controller based on the rule in Eq. (74), which can be expected in regard to the comparison of Figs. 13 and 16. As a result, the superiority of the proposed design can be summarized as follows. It leads to a decision rule on kss that guarantees rejection of disturbances due to environmental effects and inertia matrix uncertainty without unnecessarily high discontinuous control action because the complete effect of the inertia matrix uncertainty on the system is taken into account in the design. The proposed design does not necessitate trial-and-error to decide on kss . Because it uses the exact or realistic information on disturbances, which results in reasonably conservative operation of the control system, the need for utilizing adaptive sliding mode control is alleviated. © ASCE This work is supported by the Science Fellowships and Grant Programmes Department (BIDEB) of the Scientific and Technological Research Council of Turkey (TUBITAK). The authors would like to thank the dear editor and the dear reviewers for their comments and suggestions that improved the quality of the paper. References Bak, T. (1999). “Spacecraft attitude determination: A magnetometer approach.” Ph.D. thesis, Aalborg Univ., Aalborg, Denmark. CelesTrak. (2015). “NORAD two-line element sets.” 〈http://celestrak.com /NORAD/elements/〉 (May 2, 2015). Chen, Y. P., and Lo, S. C. (1993). “Sliding-mode controller design for spacecraft attitude tracking maneuvers.” IEEE Trans. Aerospace Electronic Syst., 29(4), 1328–1333. Chen, Z., Cong, B. L., and Liu, X. D. (2014). “A robust attitude control strategy with guaranteed transient performance via modified Lyapunov-based control and integral sliding mode control.” Nonlinear Dyn., 78(3), 2205–2218. Chobotov, V. A. (1991). Spacecraft attitude dynamics and control, Krieger Publishing Company, Malabar, FL. 04017082-14 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng. Downloaded from ascelibrary.org by University Of Florida on 10/25/17. Copyright ASCE. For personal use only; all rights reserved. Cong, B., Chen, Z., and Liu, X. (2012). “Disturbance observer-based adaptive integral sliding mode control for rigid spacecraft attitude maneuvers.” J. Aerosp. Eng., 227(10), 1660–1671. Davis, J. (2004). “Mathematical modeling of Earth’s magnetic field.” Technical Note, Virginia Polytechnic Institute and State Univ., Blacksburg, VA. Edwards, C., and Spurgeon, S. K. (1998). Sliding mode control: Theory and applications, Taylor & Francis, London. Fortescue, P., Swinerd, G., and Stark, J. (2011). Spacecraft systems engineering, 4th Ed., Wiley, Hoboken, NJ. Gregory, B. S. (2004). “Attitude control system design for ION, the Illinois observing nanosatellite.” M.Sc. thesis, Univ. of Illinois at UrbanaChampaign, Urbana, IL. Henderson, H. V., and Searle, S. R. (1981). “On deriving the inverse of a sum of matrices.” SIAM Rev., 23(1), 53–60. Hu, Q., Xie, L., and Gao, H. (2007). “Adaptive variable structure and active vibration reduction for flexible spacecraft under input nonlinearity.” J. Vib. Control, 13(11), 1573–1602. Hughes, P. C. (2004). Spacecraft attitude dynamics, Dover Publications, Mineola, NY. Hung, J. Y., Gao, W., and Hung, J. C. (1993). “Variable structure control: A survey.” IEEE Trans. Electron. Mag., 40(1), 2–8. Jafarov, E. M. (2009). Variable structure control and time-delay systems, WSEAS Press, Athens, Greece. Larson, W. J., and Wertz, J. R. (1999). Space mission analysis and design, 3rd Ed., Microcosm Press and Kluwer Academic Publishers, El Segundo, CA. Lo, S. C., and Chen, Y. P. (1995). “Smooth sliding-mode control for spacecraft attitude tracking maneuvers.” J. Guidance Control Dyn., 18(6), 1345–1349. Meng, Q., Zhang, T., and Song, J. Y. (2012). “Modified model-based faulttolerant time-varying attitude tracking control of uncertain flexible satellites.” J. Aerosp. Eng., 227(11), 1827–1841. © ASCE Meyer, C. D. (2000). Matrix analysis and applied linear algebra, SIAM, Philadelphia. Montenbruck, O., and Gill, E. (2000). Satellite orbits: Models, methods and applications, Springer, Berlin. Shtessel, Y., Edwards, C., Fridman, L., and Levant, A. (2014). Sliding mode control and observation, Springer, New York. Sidi, M. J. (1997). Spacecraft dynamics and control: A practical engineering approach, Cambridge University Press, New York. Sofyalı, A., Jafarov, E. M., and Wisniewski, R. (2015). “Time-varying sliding mode in rigid body motion controlled by magnetic torque.” Proc., Int. Workshop on Recent Advances in Sliding Modes (RASM 2015), Istanbul, Turkey. Utkin, V. I. (1992). Sliding modes in control and optimization, Springer, London. Vadali, S. R. (1986). “Variable-structure control of spacecraft large-angle maneuvers.” J. Guidance Control Dyn., 9(2), 235–239. Vallado, D. A. (1997). Fundamentals of astrodynamics and applications, McGraw-Hill, New York. Wertz, J. R. (1978). Spacecraft attitude determination and control, Reidel Publishing Company, Dordrecht, Netherlands. Wie, B. (1998). Space vehicle dynamics and control, American Institute of Aeronautics and Astronautics, Reston, VA. Wisniewski, R. (1996). “Satellite attitude control using only electromagnetic actuation.” Ph.D. thesis, Aalborg Univ., Aalborg, Denmark. Yang, C. D., and Sun, Y. P. (2002). “Mixed H2 =H∞ state-feedback design for microsatellite attitude control.” Control Eng. Practice, 10(9), 951–970. Zhao, J., Jiang, B., Shi, P., Gao, Z., and Xu, D. (2012). “Fault-tolerant control design for near-space vehicles based on a dynamic terminal sliding mode technique.” J. Systems and Control Eng., 226(6), 787–794. 04017082-15 J. Aerosp. Eng., 2018, 31(1): 04017082 J. Aerosp. Eng.

1/--страниц