# Super-Twisting Sliding Mode Control Law Design for Attitude Tracking Task of a Spacecraft via Reaction Wheels

Authors: Yang-Rui Li; Chao-Chung Peng
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The attitude control has been recognized as one of the most important research topics for spacecraft. If the desired attitude trajectory cannot be tracked precisely, it may cause mission failures. In the real space mission environment, the unknown external perturbations, for example, atmospheric drag and solar radiation, should be taken into consideration. Such external perturbations could deviate the precision of the spacecraft orientation and thereby lead to a mission failure. Therefore, in this paper, a quaternion-based super-twisting sliding mode robust control law for the spacecraft attitude tracking is developed. The finite time stability based on the formulation of the linear matrix inequality (LMI) is also provided. To avoid losing the control degree of freedom due to the certain actuator fault, a redundant reaction wheels configuration is adopted. The actuators distribution associated force distribution matrix (FDM) is analyzed in detail. Finally, the reference tangent-normal-binormal (TNB) command generation strategy is implemented for simulating the scenario of the space mission. Finally, the simulation results reveal that the spacecraft can achieve the desired attitude trajectory tracking demands in the presence of the time-varying external disturbances.


---

## Body

## 1. Introduction

Sliding mode control (SMC) techniques have been a popular research topic of the control theory in recent years such as adaptive super-twisting SMC [1–3], fractional-order sliding mode control [4], finite time control [5], robust backstepping SMC [6–8], and model predictive SMC [9]. The superior robustness to the matched perturbations is one of the features of the SMC. However, the price of the robustness is the chattering effect of the control signal. It causes application difficulty for practical implementations [10]. The ways to attenuate the chattering phenomenon include the following [11–13]: (i) replacing the discontinuous switching function with a saturation function or a sigmoid function, (ii) applying an adaptive law to adjust the switching gain dynamically, and (iii) using the higher-order SMC techniques. Nevertheless, skill (i) results in losing the robustness to the disturbances. Even though approach (ii) can estimate an adequate magnitude of the switching gain with respect to perturbations [14, 15], the estimation of the switch gain could increase monotonically due to the absence of perfect sliding motion in practice. For (iii), the gain/stability determinations are quite challenging.The high-order SMC approach can drive the sliding variable and its consecutive derivations to zero in the presence of the matched perturbations. However, the main challenge of the high-order SMC is that it uses the information of the high-order time derivatives of the sliding variable [16–18]. Among the higher-order SMC techniques, it is worth remarking that the second-order SMC such as the super-twisting algorithm only needs the feedback information of the sliding variable in control process. The super-twisting algorithm was firstly proposed by Dr. Levant in 1993 [19]. A quadratic Lyapunov function proposed in [20] is considered in the proof of the finite-time convergence property. The successive researches include [21–24]. Owing to the superior properties, the super-twisting algorithm has been applied in several studies, including quadrotor [25, 26], industrial emulator [27], and mobile wheeled inverted pendulum [28]. For this reason, the robust continuous super-twisting sliding mode algorithm will be adopted for the attitude tracking control design in this paper.The reaction wheel driven based system is actuated by means of generating the reaction torque from the wheel. The reaction wheel has been widely applied in the most dynamics systems. In literature [29–31], the reaction wheels are used in the attitude tracking control demands. The works in [32, 33] use the reaction wheel to address the balancing control of the inverted pendulum. The main objective of this paper is to apply the super-twisting algorithm for the attitude control of the spacecraft via using four reaction wheels as the actuated source. The attitude representation includes several approaches, eor example, Euler angles, Rodrigues parameters, and quaternion [34]. To avoid singularity, the quaternion-based control is considered. The quaternion-based control has been proposed in several studies [30, 35, 36]. However, it assumes that the scalar component of error quaternion q0e does not equal zero to guarantee that the matrix 0.5q0eI3+qe× is invertible. This assumption leads to the controllers containing a singularity when q0e=0. It should be noted that one of the reasons to use quaternion-based control is to obtain the full attitude tracking task and avoid any singularity limitations. As a result, in this paper, a quaternion-based super-twisting sliding mode algorithm is adopted in the controller design such that the robust performance can be guaranteed. The asymptotic stability proof of the nonlinear reduced-order dynamics by means of an analytic solution will be addressed without imposing assumptions.Regarding the organization of this article, in Section2, the governing equations of attitude dynamics based on the quaternion kinematics and the redundant reaction wheels configurations are derived. The configuration is introduced from [30, 37, 38]. To obtain the feasible reaction torques, the FDM is reformulated as a square and invertible matrix, which minimizes the control energy cost [30]. In Section 3, a robust, continuous super-twisting sliding mode algorithm is considered in the controller design so that the spacecraft can handle the external perturbations in the real and complex space environment. The stability problem will be reformulated as a feasibility problem of a LMI and therefore the finite time stability can be achieved in the sense of Lyapunov. In Section 4, the reference TNB command generation strategy is proposed to verify the tracking performance of the spacecraft. In Section 5, the numerical simulation is carried out and the results reveal that the spacecraft can track the desired attitude trajectory in the presence of time-varying disturbances.The contributions of this paper are summarized as follows: (i) realizes the super-twisting sliding mode algorithm as a robust, continuous quaternion-based attitude controller for the attitude trajectory tracking demands of a spacecraft with the redundant reaction wheels; (ii) proposes a modified version of LMI which has higher degrees of freedom for finding the decision variables and it can satisfy the convergence performance by requirement; (iii) derives the analytic solution of the nonlinear reduced-order dynamics; and (iv) presents a reference TNB command generation strategy so that the feasibility of the controller can be verified.


## 2. Reaction Wheels Driven Based on Spacecraft Dynamics Modeling

### 2.1. Geometry Configuration Analysis

From the perspective of practical realization, to avoid losing a degree of freedom of control in space due to certain actuator faults, the redundant reaction wheels configuration is adopted [30, 37, 38]. The dynamic configuration is shown in Figure 1. To formally derive the governing equations of the attitude dynamics, we firstly define the coordinate system as follows: (i) the body frame denoted as xyz, which is fixed in the body of the spacecraft to represent the attitude of the spacecraft, and (ii) the auxiliary rotation frame denoted as xiyizi, which is fixed in the i-th reaction wheel to describe the relative rotation of i-th reaction wheel to spacecraft.Figure 1
Geometry prototype of the redundant reaction wheels configuration.Taking the first reaction wheel as an example, the geometry mapping relation betweenxyz and x1y1z1 is explained in Figures 2 and 3, respectively. Referring to Figure 2, the body frame xyz rotates about the z-axis with an angle β, and the new frame is denoted as x′y′z′. From the rotation property, the mapping relation between x′y′z′ and xyz can be constructed as(1)x′y′z′=cosβsinβ0−sinβcosβ0001xyz.Figure 2
Rotation of body framexyz to frame x′y′z′.Figure 3
Rotation of framex′y′z′ to frame x″y″z″.Regarding Figure3, the frame x′y′z′ rotates about negative y-axis with an angle α, and then the new frame x′′y″z′′ is obtained. Again, from the rotation property, we have the mapping relation between x′y′z′ and x′′y″z′′:(2)x′′y′′z′′=cos−α0−sin−α010sin−α0cos−αx′y′z′.On the basis of previous illustrations, the mapping relation betweenxyz and x1y1z1 can be derived by combining (1) and (2):(3)x1y1z1=x′′y′′z′′=cos−α0−sin−α010sin−α0cos−αx′y′z′=cos−α0−sin−α010sin−α0cos−αcosβsinβ0−sinβcosβ0001xyz.In general, the mapping relation betweenxyz and xiyizi can be formulated as(4)xiyizi=Rixyz,where(5)Ri=cos−α0−sin−α010sin−α0cos−αcosβisinβi0−sinβicosβi0001=cosαcosβicosαsinβisinα−sinβicosβi0−cosβisinα−sinαsinβicosα,and βi=β+0.5i−1π, i=1,2,3,4.


### 2.2. Attitude Dynamics

To avoid the singularity problem in the Euler approaches, the quaternion-based attitude representation is considered. LetQ=q0,qTT=q0,q1,q2,q3T∈ℝ4 be the unit quaternion and let ω=ωx,ωy,ωzT∈ℝ3 be the angular velocity of spacecraft; then the quaternion kinematic equations are given by(6)q˙0=−12qTω,(7)q.=12q0I3+q×ω.The above equations can be rewritten in a compact matrix form:(8)Q.=12EQω,where(9)EQ=−qTq0I3+q×=−q1−q2−q3q0−q3q2q3q0−q1−q2q1q0.The symbola× represents a cross product matrix of vector a=a1,a2,a3T∈ℝ3, which is defined as(10)a×=0−a3a2a30−a1−a2a10.In the following, the kinetic equations based on the redundant reaction wheels configuration will be derived. LetΩi, i=1,2,3,4, be the relative rotational speed of each reaction wheel to the spacecraft; then the angular velocity of each reaction wheel with respect to xiyizi is given by(11)ωw,i=Riω+Ωiexi,exi=100T.The angular momentum of spacecraft relative the body framexyz is(12)H=Jω=Jx−Jxy−Jxz−JxyJy−Jyz−Jxz−JyzJzωxωyωz.Assume that the reaction wheels are perfect circle plate. Hence, the moment of inertia matrix ofi-th reaction wheel with respect to xiyizi is(13)Jw,i=Jm,i000Ji′000Ji′.The angular momentum of each reaction wheel relative framexiyizi can be constructed:(14)Hw,i=Jw,iωw,i=Jw,iRiω+Ωiexi=Jw,iRiω+JwΩiexi=Jw,iRiω+JmiΩiexi.MappingHw,i onto spacecraft body frame xyz yields(15)Hw,iB=RiTωw,i=RiTJw,iRiω+JmiΩiexi=RiTJw,iRiω+RiTJmiΩiexi.The total angular momentum of the system (spacecraft and reaction wheels),HT, is the summation of H and Hw,iB, which is(16)HT=H+∑i=14Hw,iB=Jω+∑i=14RiTJw,iRiω+RiTJmiΩiexi=J+∑i=14RiTJw,iRiω+∑i=14RiTJmiΩiexi.The term(17)∑i=14RiTJmΩiexi=cosαcosβcosαsinβsinαJm1Ω1+−cosαsinβcosαcosβsinαJm2Ω2+−cosαcosβ−cosαsinβsinαJm3Ω3+cosαcosβ−cosαsinβsinαJm4Ω4=cosαcosβ−cosαsinβ−cosαcosβcosαsinβcosαsinβcosαcosβ−cosαsinβ−cosαcosβsinαsinαsinαsinαJm10000Jm20000Jm30000Jm4Ω1Ω2Ω3Ω4.Define the following:(i)
Equivalent moment of inertia matrix of the system:(18)Jeq=J+∑i=14RiTJw,iRi.(ii)
Force distribution matrix (FDM):(19)Γ=cosαcosβ−cosαsinβ−cosαcosβcosαsinβcosαsinβcosαcosβ−cosαsinβ−cosαcosβsinαsinαsinαsinα.(iii)
Axial moment of inertia matrix:(20)Jm=Jm10000Jm20000Jm30000Jm4.Combining (17)–(20), equation (16) can be further simplified to(21)HT=Jeqω+ΓJmΩ,where Ω=Ω1,Ω2,Ω3,Ω4T. The Euler equation of motion is given by(22)MG=dHTdtB+ω×HT,in which the total external torque MG is equal to the summation of the external control torques τa and the external disturbance torques dt. Substituting (21) into (22) yields(23)τa+d=Jeqω.+ΓJmΩ.+ω×Jeqω+ΓJmΩ.Define the control torqueτ and the reaction torque τw as(24)τc=Γτw≜τcxτcyτczT,τw=−JmΩ.=−Jm1Ω˙1Jm2Ω˙2Jm3Ω˙3Jm4Ω˙4T≜τw1τw2τw3τw4T.Hence, (23) can be further simplified as(25)Jeqω.=τa+τc+d−ω×Jeqω+ΓJmΩ.In this paper,τa=0 is considered.


### 2.3. Actuator Analysis

From (24), we have(26)τcxτcyτcz=cosαcosβ−cosαsinβ−cosαcosβcosαsinβcosαsinβcosαcosβ−cosαsinβ−cosαcosβsinαsinαsinαsinατw1τw2τw3τw4.Onceτc is designed for the attitude trajectory tracking demands of the system dynamics (6), (7), and (25), it is desired to obtain each reaction torque τwi. However, FDM (19) is a nonsquare matrix, so the inverse does not exist. To obtain the FDM with a special matrix structure, the two following geometry constraints are imposed [30]:(27)sinβ=cosβ,0≤β<π2,cosαsinβ=sinα,0≤α<π2,or, equivalently,(28)β=π4,sinα=33.From (27), we have the following FDM:(29)Γ=331−1−1111−1−11111,and then the following static optimization problem is formulated [30]:(30)minτwT=∑i=14τwi2,subject to(31)g1≜33τw1−τw2−τw3+τw4−τcx=0,g2≜33τw1+τw2−τw3−τw4−τcy=0,g3≜33τw1+τw2+τw3+τw4−τcz=0.To formally address the problem, refer to [39], and define the Lagrangian ℒ together with the Lagrange multiplier λ=λ1,λ2,λ3T∈ℝ3 as follows:(32)ℒ=T+λ1g1+λ2g2+λ3g3.Letτw∗ and λ∗ be the optimal solution. The first-order necessary condition to minimize T is(33)∂L∂τw1|τ∗,λ∗=2τw1∗+33λ1∗+33λ2∗+33λ3∗=0,∂L∂τw2|τ∗,λ∗=2τw2∗−33λ1∗+33λ2∗+33λ3∗=0,∂L∂τw3|τ∗,λ∗=2τw3∗−33λ1∗−33λ2∗+33λ3∗=0,∂L∂τw4|τ∗,λ∗=2τw4∗+33λ1∗−33λ2∗+33λ3∗=0,and(34)∂L∂λ1|τ∗,λ∗=33τw1∗−τw2∗−τw3∗+τw4∗−τcx=0,∂L∂λ2|τ∗,λ∗=33τw1∗+τw2∗−τw3∗−τw4∗−τcy=0,∂L∂λ3|τ∗,λ∗=33τw1∗+τw2∗+τw3∗+τw4∗−τcz=0.From (33), it is implied that(35)τw1∗−τw2∗+τw3∗−τw4∗=0.Combining (34) and (35), the distribution matrix can be augmented as(36)τcxτcyτcz0=33−33−33333333−33−33333333331−11−1τw1∗τw2∗τw3∗τw4∗.The inverse mapping is(37)τw1∗τw2∗τw3∗τw4∗=143331−333−1−3−3313−33−1τcxτcyτcz0.Hence, whenτc is designed for the system dynamics (6), (7), and (25), the reaction torque τwi for each reaction wheel can be obtained by inverse mapping (37). Moreover, it can be guaranteed that the reaction torque is an optimal value τw∗ to minimize the performance index (30).Because the reaction wheels are actuated by servo motors, that is, the reaction torques are generated by servo motors, the following linear dynamic equation is considered:(38)JmiΩ˙i=−aiΩi+biVin,i,i=1,2,3,4.From the definition, (24) and (38) become(39)Vin,i=1bi−τwi+aiΩii=1,2,3,4,where the parameters ai=0.02 and bi=3 are considered. Vin,i is the control input voltage.Conclusively,τc is first designed for the system dynamics (6), (7), and (25). Secondly, inverse mapping (37) is used to obtain the optimal reaction torque τwi of each reaction wheel to minimize the energy cost (30). Finally, apply (39) to obtain the corresponding control voltage Vin,i. In the following, the discussion focuses on how to design τc so that the desired trajectory can be achieved in the presence of time-varying disturbances.


## 2.1. Geometry Configuration Analysis

From the perspective of practical realization, to avoid losing a degree of freedom of control in space due to certain actuator faults, the redundant reaction wheels configuration is adopted [30, 37, 38]. The dynamic configuration is shown in Figure 1. To formally derive the governing equations of the attitude dynamics, we firstly define the coordinate system as follows: (i) the body frame denoted as xyz, which is fixed in the body of the spacecraft to represent the attitude of the spacecraft, and (ii) the auxiliary rotation frame denoted as xiyizi, which is fixed in the i-th reaction wheel to describe the relative rotation of i-th reaction wheel to spacecraft.Figure 1
Geometry prototype of the redundant reaction wheels configuration.Taking the first reaction wheel as an example, the geometry mapping relation betweenxyz and x1y1z1 is explained in Figures 2 and 3, respectively. Referring to Figure 2, the body frame xyz rotates about the z-axis with an angle β, and the new frame is denoted as x′y′z′. From the rotation property, the mapping relation between x′y′z′ and xyz can be constructed as(1)x′y′z′=cosβsinβ0−sinβcosβ0001xyz.Figure 2
Rotation of body framexyz to frame x′y′z′.Figure 3
Rotation of framex′y′z′ to frame x″y″z″.Regarding Figure3, the frame x′y′z′ rotates about negative y-axis with an angle α, and then the new frame x′′y″z′′ is obtained. Again, from the rotation property, we have the mapping relation between x′y′z′ and x′′y″z′′:(2)x′′y′′z′′=cos−α0−sin−α010sin−α0cos−αx′y′z′.On the basis of previous illustrations, the mapping relation betweenxyz and x1y1z1 can be derived by combining (1) and (2):(3)x1y1z1=x′′y′′z′′=cos−α0−sin−α010sin−α0cos−αx′y′z′=cos−α0−sin−α010sin−α0cos−αcosβsinβ0−sinβcosβ0001xyz.In general, the mapping relation betweenxyz and xiyizi can be formulated as(4)xiyizi=Rixyz,where(5)Ri=cos−α0−sin−α010sin−α0cos−αcosβisinβi0−sinβicosβi0001=cosαcosβicosαsinβisinα−sinβicosβi0−cosβisinα−sinαsinβicosα,and βi=β+0.5i−1π, i=1,2,3,4.


## 2.2. Attitude Dynamics

To avoid the singularity problem in the Euler approaches, the quaternion-based attitude representation is considered. LetQ=q0,qTT=q0,q1,q2,q3T∈ℝ4 be the unit quaternion and let ω=ωx,ωy,ωzT∈ℝ3 be the angular velocity of spacecraft; then the quaternion kinematic equations are given by(6)q˙0=−12qTω,(7)q.=12q0I3+q×ω.The above equations can be rewritten in a compact matrix form:(8)Q.=12EQω,where(9)EQ=−qTq0I3+q×=−q1−q2−q3q0−q3q2q3q0−q1−q2q1q0.The symbola× represents a cross product matrix of vector a=a1,a2,a3T∈ℝ3, which is defined as(10)a×=0−a3a2a30−a1−a2a10.In the following, the kinetic equations based on the redundant reaction wheels configuration will be derived. LetΩi, i=1,2,3,4, be the relative rotational speed of each reaction wheel to the spacecraft; then the angular velocity of each reaction wheel with respect to xiyizi is given by(11)ωw,i=Riω+Ωiexi,exi=100T.The angular momentum of spacecraft relative the body framexyz is(12)H=Jω=Jx−Jxy−Jxz−JxyJy−Jyz−Jxz−JyzJzωxωyωz.Assume that the reaction wheels are perfect circle plate. Hence, the moment of inertia matrix ofi-th reaction wheel with respect to xiyizi is(13)Jw,i=Jm,i000Ji′000Ji′.The angular momentum of each reaction wheel relative framexiyizi can be constructed:(14)Hw,i=Jw,iωw,i=Jw,iRiω+Ωiexi=Jw,iRiω+JwΩiexi=Jw,iRiω+JmiΩiexi.MappingHw,i onto spacecraft body frame xyz yields(15)Hw,iB=RiTωw,i=RiTJw,iRiω+JmiΩiexi=RiTJw,iRiω+RiTJmiΩiexi.The total angular momentum of the system (spacecraft and reaction wheels),HT, is the summation of H and Hw,iB, which is(16)HT=H+∑i=14Hw,iB=Jω+∑i=14RiTJw,iRiω+RiTJmiΩiexi=J+∑i=14RiTJw,iRiω+∑i=14RiTJmiΩiexi.The term(17)∑i=14RiTJmΩiexi=cosαcosβcosαsinβsinαJm1Ω1+−cosαsinβcosαcosβsinαJm2Ω2+−cosαcosβ−cosαsinβsinαJm3Ω3+cosαcosβ−cosαsinβsinαJm4Ω4=cosαcosβ−cosαsinβ−cosαcosβcosαsinβcosαsinβcosαcosβ−cosαsinβ−cosαcosβsinαsinαsinαsinαJm10000Jm20000Jm30000Jm4Ω1Ω2Ω3Ω4.Define the following:(i)
Equivalent moment of inertia matrix of the system:(18)Jeq=J+∑i=14RiTJw,iRi.(ii)
Force distribution matrix (FDM):(19)Γ=cosαcosβ−cosαsinβ−cosαcosβcosαsinβcosαsinβcosαcosβ−cosαsinβ−cosαcosβsinαsinαsinαsinα.(iii)
Axial moment of inertia matrix:(20)Jm=Jm10000Jm20000Jm30000Jm4.Combining (17)–(20), equation (16) can be further simplified to(21)HT=Jeqω+ΓJmΩ,where Ω=Ω1,Ω2,Ω3,Ω4T. The Euler equation of motion is given by(22)MG=dHTdtB+ω×HT,in which the total external torque MG is equal to the summation of the external control torques τa and the external disturbance torques dt. Substituting (21) into (22) yields(23)τa+d=Jeqω.+ΓJmΩ.+ω×Jeqω+ΓJmΩ.Define the control torqueτ and the reaction torque τw as(24)τc=Γτw≜τcxτcyτczT,τw=−JmΩ.=−Jm1Ω˙1Jm2Ω˙2Jm3Ω˙3Jm4Ω˙4T≜τw1τw2τw3τw4T.Hence, (23) can be further simplified as(25)Jeqω.=τa+τc+d−ω×Jeqω+ΓJmΩ.In this paper,τa=0 is considered.


## 2.3. Actuator Analysis

From (24), we have(26)τcxτcyτcz=cosαcosβ−cosαsinβ−cosαcosβcosαsinβcosαsinβcosαcosβ−cosαsinβ−cosαcosβsinαsinαsinαsinατw1τw2τw3τw4.Onceτc is designed for the attitude trajectory tracking demands of the system dynamics (6), (7), and (25), it is desired to obtain each reaction torque τwi. However, FDM (19) is a nonsquare matrix, so the inverse does not exist. To obtain the FDM with a special matrix structure, the two following geometry constraints are imposed [30]:(27)sinβ=cosβ,0≤β<π2,cosαsinβ=sinα,0≤α<π2,or, equivalently,(28)β=π4,sinα=33.From (27), we have the following FDM:(29)Γ=331−1−1111−1−11111,and then the following static optimization problem is formulated [30]:(30)minτwT=∑i=14τwi2,subject to(31)g1≜33τw1−τw2−τw3+τw4−τcx=0,g2≜33τw1+τw2−τw3−τw4−τcy=0,g3≜33τw1+τw2+τw3+τw4−τcz=0.To formally address the problem, refer to [39], and define the Lagrangian ℒ together with the Lagrange multiplier λ=λ1,λ2,λ3T∈ℝ3 as follows:(32)ℒ=T+λ1g1+λ2g2+λ3g3.Letτw∗ and λ∗ be the optimal solution. The first-order necessary condition to minimize T is(33)∂L∂τw1|τ∗,λ∗=2τw1∗+33λ1∗+33λ2∗+33λ3∗=0,∂L∂τw2|τ∗,λ∗=2τw2∗−33λ1∗+33λ2∗+33λ3∗=0,∂L∂τw3|τ∗,λ∗=2τw3∗−33λ1∗−33λ2∗+33λ3∗=0,∂L∂τw4|τ∗,λ∗=2τw4∗+33λ1∗−33λ2∗+33λ3∗=0,and(34)∂L∂λ1|τ∗,λ∗=33τw1∗−τw2∗−τw3∗+τw4∗−τcx=0,∂L∂λ2|τ∗,λ∗=33τw1∗+τw2∗−τw3∗−τw4∗−τcy=0,∂L∂λ3|τ∗,λ∗=33τw1∗+τw2∗+τw3∗+τw4∗−τcz=0.From (33), it is implied that(35)τw1∗−τw2∗+τw3∗−τw4∗=0.Combining (34) and (35), the distribution matrix can be augmented as(36)τcxτcyτcz0=33−33−33333333−33−33333333331−11−1τw1∗τw2∗τw3∗τw4∗.The inverse mapping is(37)τw1∗τw2∗τw3∗τw4∗=143331−333−1−3−3313−33−1τcxτcyτcz0.Hence, whenτc is designed for the system dynamics (6), (7), and (25), the reaction torque τwi for each reaction wheel can be obtained by inverse mapping (37). Moreover, it can be guaranteed that the reaction torque is an optimal value τw∗ to minimize the performance index (30).Because the reaction wheels are actuated by servo motors, that is, the reaction torques are generated by servo motors, the following linear dynamic equation is considered:(38)JmiΩ˙i=−aiΩi+biVin,i,i=1,2,3,4.From the definition, (24) and (38) become(39)Vin,i=1bi−τwi+aiΩii=1,2,3,4,where the parameters ai=0.02 and bi=3 are considered. Vin,i is the control input voltage.Conclusively,τc is first designed for the system dynamics (6), (7), and (25). Secondly, inverse mapping (37) is used to obtain the optimal reaction torque τwi of each reaction wheel to minimize the energy cost (30). Finally, apply (39) to obtain the corresponding control voltage Vin,i. In the following, the discussion focuses on how to design τc so that the desired trajectory can be achieved in the presence of time-varying disturbances.


## 3. Super-Twisting Sliding Mode Controller Design

### 3.1. Super-Twisting Sliding Mode Algorithm

The design process of sliding mode control includes two steps: (i) A sliding variables is designed so that the stability of the reduced-order dynamics can be guaranteed. (ii) Seeking the robust, continuous control law to guarantee the sliding mode s=s˙=0 occurs in a finite time in the presence of the time-varying external disturbances. In the following, the super-twisting sliding mode algorithm [19] is introduced.The super-twisting sliding mode algorithm is one of the outstanding robust control algorithms which handles a system with a relative degree (u⟶s) equal to one. Based on the algorithm, the closed-loop sliding dynamics is designed as(40)s˙=−k1s1/2signs−k2∫0tsignsdτ+dt,where s=sx,t∈ℝ is the sliding variable; x∈ℝn is state vector; dt∈ℝ is an unknown time-varying perturbation and it is assumed that d˙t≤δ; the gain pair k1,k2 is to be designed so that the sliding mode can occur in a finite time.In fact, the more general representation of (40) is to express it as the following state-space form. Let(41)z1=sx,t.z2=−k2∫0tsignsdτ+dt,which implies(42)z˙1=−k1z11/2signz1+z2,z˙2=−k2signz1+ρ,where ρt=d˙t and ρt=d˙t≤δ. Since (42) is nonlinear, consider the following variable transformation [20]:(43)ζ1=z11/2signz1,ζ2=z2.Taking the time derivative yields(44)ζ˙1=z˙1z11/2−z1⋅1/2z1−1/2⋅signz1⋅z˙1z1=−k1z11/2signz1+z2⋅z11/2−1/2z11/2z1=1z11/2−k12z11/2signz1+12z2=1ζ1−k12ζ1+12ζ2,ζ˙2=−k2signz1+ρ=1z11/2−k2z11/2signz1+z11/2ρ=1ζ1−k2ζ1+ζ1ρ,which can be rewritten in the matrix form as(45)ζ.=1ζ1Aζ+Bρ˜,A=−0.5k10.5−k20,B=01,where ζ=ζ1,ζ2T and ζ1=z11/2. Disturbance transformation ρ˜t,ζ1=ρtζ1 satisfies(46)ρ˜t,ζ1≤δζ1.For system (45), the stability problem is proven in the following. It is shown that the stability issue can be reformulated as a feasibility problem in terms of the LMI.Theorem 1.
(See[21, 22, 25]). Suppose that there exist symmetric and positive definite matrices P=PT>0 and Qc=QcT>0 so that the following LMI,(47)PA+ATP+Qc+CTCPBBTP−γ2<0,is feasible, where γ=1/δ, C=10, and A and B are provided in (45). Then the quadratic form,(48)V=ζTPζ,is a strict Lyapunov function for system (45) and the trajectory reaches the origin in a finite time.Proof of Theorem 1.
Applying the Rayleigh inequality,V is bounded by(49)λminPζ2≤V≤λmaxPζ2,where ζ2=z1+z22 represents the Euclidean norm of ζ. For (48), taking time derivative gives(50)V˙=1ζ1ζTPA+ATPζ+ζTPBρ˜+ρ˜BTPζ.
According to (46), the following inequality is satisfied:(51)ρ˜t,ζ12≤δ2ζ12≤δ2ζ12+ζ22,which guarantees(52)ζTζ−1δ2ρ˜t,ζ12>0.
Hence, (50) can be rewritten as(53)V˙≤1ζ1ζTPA+ATPζ+ζTPBρ˜+ρ˜BTPζ+ζTζ−1δ2ρ˜2=1ζ1ζTPA+ATP+CTC+Qc−Qcζ+ζTPBρ˜+ρ˜BTPζ−1δ2ρ˜2=−1ζ1ζTQcζ+1ζ1ζρ˜TPA+ATP+Qc+CTCPBBTP−γ2ζρ˜.where γ=1/δ.
Based on (47) and (49), (53) can be further simplified to(54)V˙≤−1ζ1ζTQcζ≤−1ζ1λminQcζ2≤−1ζ1λminQcVλmaxP.
Again, from (49), the following inequality can be deduced:(55)−1ζ1≤−1ζ≤−λmin1/2PV1/2.
Based on (55), it can be concluded that (54) satisfies(56)V˙≤−αV1/2,α=λminQcλmin1/2PλmaxP.
For (56), one has(57)∫V0VtdVV1/2≤−α∫0tdt,which implies(58)Vt≤V1/20−α2t2.
Hence,Vt reaches zero within a finite time described by(59)tf≤2V1/20α,where V0 is the initial condition of Vt and α is given by (56). Based on the appropriate gain pair k1,k2 selection, the LMI equation (47) can be established. As a result, the finite time stability as shown by (59) can be achieved.
In Theorem1, we suppose that the LMI (47) is negative definite. Theorem 2 is similar to [21, 22, 25] and is presented for the feasibility of LMI (48) in Theorem 1.Theorem 2.
Consider the LMI given by (47); there exists a feasible solution P, Qc so that the LMI (47) can be established if and only if the parameters k1 and k2 in A satisfy(60)k2>δ,k12>4k2,or (61)k1212k2−116k12<δ2,4k2>k12.
Furthermore, the additional constraintsk2≠0 for (60) and k1≠0 and k12≠8k2 for (61) must also be satisfied.Proof of Theorem 2.
If the LMI (45) is feasible, then the ℒ2-gain of the following system,(62)Gs=1/2s2+1/2k1s+1/2k2,must be less than or equal to γ; that is,(63)maxωGjω<γ=1δ⇒maxωGjω2<1δ2.
The above statement is the so-called bounded-real condition [40]. In order to find proper k1,k2 in (62) so that condition (63) can be satisfied, calculate(64)Gjω2=1k2−2ω22+k1ω2,and its derivative(65)ddωGjω2=−16ωω2+1/8k12−1/2k2k2−2ω22+k1ω22.
The extreme point can be obtained by setting (65) equal to zero. Checking for the second-order sufficient condition d2/dω2Gjω2, it can be deduced that maxωGjω can be reached, when(66)ω=0,if4k2−k12<0,4k2−k1281/2,if4k2−k12>0.
Substituting (66) into (64) yields(67)maxωGjω2=1k22,if4k2−k12<0.1k121/2k2−1/16k12,if4k2−k12>0.
Thus, combining (63) and (67) shows that if(68)k2>δ,k12>4k2,(69)ork1212k2−116k12>δ2,4k2>k12,then the LMI (47) is feasible. Moreover, to avoid the singularity of (67), the constraints k2≠0 for (68) and k1≠0 and k12≠8k2 for (69) are made.Remark 1.
Notice that system (62) is not a transfer function of system (45); it is the corresponding linear system of LMI (47). More details can be found in [40].


### 3.2. Controller Design

The control objective is to design a robust control torqueτc such that the spacecraft can achieve the arbitrary attitude trajectory tracking demands in the presence of time-varying disturbances.For this reason, letQd=q0d,qdTT=q0d,q1d,q2d,q3dT∈ℝ4 be the desired quaternion and let ωd=ωxd,ωyd,ωzdT∈ℝ3 be the desired angular velocity. Define the tracking error vectors Qe=q0e,qeTT∈ℝ4 and ωe=ωxe,ωye,ωzeT∈ℝ3 as follows [36]:(70)Qe=Qd−1⊗Q=q0eqe=q0dq0+qdTqq0dq−q0qd−qd×q,ωe=ω−ωd=ωxeωyeωze=ωx−ωxdωy−ωydωz−ωzd,where ⊗ represents the quaternion multiplication. Based on (70), the quaternion-based error dynamics of (6), (7), and (25) is given by(71)q˙0e=−12qeTωe,(72)q˙e=12q0eI3+qe×ωe,(73)ω.e=Jeq−1−τc+d−ω×Jeqω+ΓJmΩ−ω.d.Select the sliding surface as(74)S=ωe+λqe,where S=s1,s2,s3T∈ℝ3; λ∈ℝ1 is a positive parameter to be designed.Suppose that the sliding motion is fulfilled in a finite timet=tf; it gives(75)S=S˙=0,∀t≥tf.From (71), (72), (74), and (75), the nonlinear reduced-order dynamics can be obtained:(76)q˙0e=−12qeTωe=−12qeT−λqe=12λqeTqe,∀t≥tf,(77)q˙e=12q0eI3+qe×−λqe=−12λq0eqe,∀t≥tf.The term “reduced-order” means that the system error dynamics described by (71)–(73 with order 7 “reduced” to subsystem (76) and (77) with order 4. In order to analyze the stability of the reduced-order dynamics, apply the identity of unit quaternion:(78)q0e2+qeTqe=1.Equation (76) can be decoupled as(79)q˙0e=12λ1−q0e2,∀t≥tf.Then, consider the following variable transformation:(80)q0et=1+1yt,which implies(81)q˙0et=−y˙ty2.Substituting (79) yields(82)−y˙ty2=12λ1−1+1yt2=12λ−2y−1y2,which implies(83)y˙−λy=12λ.Clearly, equation (83) is a linear ordinary differential equation. Its solution is given by(84)yt=Ceλt−12,where C is an integration constant. Apply the inverse mapping from (80); we get(85)q0et=1+1−1/2+Ceλt.Letq0e0 be the initial condition; we have(86)C=12+1q0e0−1.Hence, the analytic solution of (79) can be obtained:(87)q0et=1+2q0etf−1−q0etf+1+q0etf+1eλt−tf,∀t≥tf.Observing (87), it can be deduced that (i) the singularity occurs as an improper sliding gain λ is chosen such that −q0etf+1+q0etf+1eλt−tf=0,∀t≥tf. Thus, the sliding gain λ>0 is designed so that −q0etf+1+q0etf+1eλt−tf≠0,∀t≥tf. (ii) q0e⟶+1 as t⟶∞ from (87) and qe⟶0 as t⟶∞ from the identity of the unit quaternion (78). Observing (70), it can be found that q0e⟶1 and qe⟶0 imply q0⟶q0d and q⟶⟶q⟶d. That is, the nonlinear reduced-order dynamics (76) and (77) are asymptotically stable. It is different and more outstanding than [30, 35, 36]; the assumption that q0e≠0 is not made. In the following, how to enter the sliding mode in a finite time in the presence of external disturbances by means of the robust control law will be discussed.In order to introduceτc, taking the time derivative about (74) yields(88)S˙=ω.e+λq˙e=Jeq−1−τc+d−ω×Jeqω+ΓJmΩ+ξ,where(89)ξ=−ω.d+0.5λq0eI3+qe×ωe.Based on the super-twisting sliding mode algorithm [21], the following robust control law is designed:(90)τc0=ω×Jeqω+ΓJmΩ−Jeqξ,τcN=−JeqK1SS1/2+K2∫0tsignSτdτ,τc=−τc0+τcN,in which the gain matrices are denoted as(91)K1=diagk11k12k13,K2=diagk21k22k23.Substituting (90) into (88) yields the closed-loop sliding dynamics:(92)S˙=−K1SS1/2−K2∫0tsignSτdτ+D,where D=Jeq−1d=D1,D2,D3T∈ℝ3 and signS∈ℝ3 is a sign function defined as follows:(93)signsi=1,ifsi>0,−1,ifsi<0,signsi∈−1,1,ifsi=0,i=1,2,3..According to the stability criteria derived from Theorems1 and 2, the following gains are chosen:(94)k2i>δi,k1i>4k2i,where δi=supD˙i, i=1,2,3.


## 3.1. Super-Twisting Sliding Mode Algorithm

The design process of sliding mode control includes two steps: (i) A sliding variables is designed so that the stability of the reduced-order dynamics can be guaranteed. (ii) Seeking the robust, continuous control law to guarantee the sliding mode s=s˙=0 occurs in a finite time in the presence of the time-varying external disturbances. In the following, the super-twisting sliding mode algorithm [19] is introduced.The super-twisting sliding mode algorithm is one of the outstanding robust control algorithms which handles a system with a relative degree (u⟶s) equal to one. Based on the algorithm, the closed-loop sliding dynamics is designed as(40)s˙=−k1s1/2signs−k2∫0tsignsdτ+dt,where s=sx,t∈ℝ is the sliding variable; x∈ℝn is state vector; dt∈ℝ is an unknown time-varying perturbation and it is assumed that d˙t≤δ; the gain pair k1,k2 is to be designed so that the sliding mode can occur in a finite time.In fact, the more general representation of (40) is to express it as the following state-space form. Let(41)z1=sx,t.z2=−k2∫0tsignsdτ+dt,which implies(42)z˙1=−k1z11/2signz1+z2,z˙2=−k2signz1+ρ,where ρt=d˙t and ρt=d˙t≤δ. Since (42) is nonlinear, consider the following variable transformation [20]:(43)ζ1=z11/2signz1,ζ2=z2.Taking the time derivative yields(44)ζ˙1=z˙1z11/2−z1⋅1/2z1−1/2⋅signz1⋅z˙1z1=−k1z11/2signz1+z2⋅z11/2−1/2z11/2z1=1z11/2−k12z11/2signz1+12z2=1ζ1−k12ζ1+12ζ2,ζ˙2=−k2signz1+ρ=1z11/2−k2z11/2signz1+z11/2ρ=1ζ1−k2ζ1+ζ1ρ,which can be rewritten in the matrix form as(45)ζ.=1ζ1Aζ+Bρ˜,A=−0.5k10.5−k20,B=01,where ζ=ζ1,ζ2T and ζ1=z11/2. Disturbance transformation ρ˜t,ζ1=ρtζ1 satisfies(46)ρ˜t,ζ1≤δζ1.For system (45), the stability problem is proven in the following. It is shown that the stability issue can be reformulated as a feasibility problem in terms of the LMI.Theorem 1.
(See[21, 22, 25]). Suppose that there exist symmetric and positive definite matrices P=PT>0 and Qc=QcT>0 so that the following LMI,(47)PA+ATP+Qc+CTCPBBTP−γ2<0,is feasible, where γ=1/δ, C=10, and A and B are provided in (45). Then the quadratic form,(48)V=ζTPζ,is a strict Lyapunov function for system (45) and the trajectory reaches the origin in a finite time.Proof of Theorem 1.
Applying the Rayleigh inequality,V is bounded by(49)λminPζ2≤V≤λmaxPζ2,where ζ2=z1+z22 represents the Euclidean norm of ζ. For (48), taking time derivative gives(50)V˙=1ζ1ζTPA+ATPζ+ζTPBρ˜+ρ˜BTPζ.
According to (46), the following inequality is satisfied:(51)ρ˜t,ζ12≤δ2ζ12≤δ2ζ12+ζ22,which guarantees(52)ζTζ−1δ2ρ˜t,ζ12>0.
Hence, (50) can be rewritten as(53)V˙≤1ζ1ζTPA+ATPζ+ζTPBρ˜+ρ˜BTPζ+ζTζ−1δ2ρ˜2=1ζ1ζTPA+ATP+CTC+Qc−Qcζ+ζTPBρ˜+ρ˜BTPζ−1δ2ρ˜2=−1ζ1ζTQcζ+1ζ1ζρ˜TPA+ATP+Qc+CTCPBBTP−γ2ζρ˜.where γ=1/δ.
Based on (47) and (49), (53) can be further simplified to(54)V˙≤−1ζ1ζTQcζ≤−1ζ1λminQcζ2≤−1ζ1λminQcVλmaxP.
Again, from (49), the following inequality can be deduced:(55)−1ζ1≤−1ζ≤−λmin1/2PV1/2.
Based on (55), it can be concluded that (54) satisfies(56)V˙≤−αV1/2,α=λminQcλmin1/2PλmaxP.
For (56), one has(57)∫V0VtdVV1/2≤−α∫0tdt,which implies(58)Vt≤V1/20−α2t2.
Hence,Vt reaches zero within a finite time described by(59)tf≤2V1/20α,where V0 is the initial condition of Vt and α is given by (56). Based on the appropriate gain pair k1,k2 selection, the LMI equation (47) can be established. As a result, the finite time stability as shown by (59) can be achieved.
In Theorem1, we suppose that the LMI (47) is negative definite. Theorem 2 is similar to [21, 22, 25] and is presented for the feasibility of LMI (48) in Theorem 1.Theorem 2.
Consider the LMI given by (47); there exists a feasible solution P, Qc so that the LMI (47) can be established if and only if the parameters k1 and k2 in A satisfy(60)k2>δ,k12>4k2,or (61)k1212k2−116k12<δ2,4k2>k12.
Furthermore, the additional constraintsk2≠0 for (60) and k1≠0 and k12≠8k2 for (61) must also be satisfied.Proof of Theorem 2.
If the LMI (45) is feasible, then the ℒ2-gain of the following system,(62)Gs=1/2s2+1/2k1s+1/2k2,must be less than or equal to γ; that is,(63)maxωGjω<γ=1δ⇒maxωGjω2<1δ2.
The above statement is the so-called bounded-real condition [40]. In order to find proper k1,k2 in (62) so that condition (63) can be satisfied, calculate(64)Gjω2=1k2−2ω22+k1ω2,and its derivative(65)ddωGjω2=−16ωω2+1/8k12−1/2k2k2−2ω22+k1ω22.
The extreme point can be obtained by setting (65) equal to zero. Checking for the second-order sufficient condition d2/dω2Gjω2, it can be deduced that maxωGjω can be reached, when(66)ω=0,if4k2−k12<0,4k2−k1281/2,if4k2−k12>0.
Substituting (66) into (64) yields(67)maxωGjω2=1k22,if4k2−k12<0.1k121/2k2−1/16k12,if4k2−k12>0.
Thus, combining (63) and (67) shows that if(68)k2>δ,k12>4k2,(69)ork1212k2−116k12>δ2,4k2>k12,then the LMI (47) is feasible. Moreover, to avoid the singularity of (67), the constraints k2≠0 for (68) and k1≠0 and k12≠8k2 for (69) are made.Remark 1.
Notice that system (62) is not a transfer function of system (45); it is the corresponding linear system of LMI (47). More details can be found in [40].


## 3.2. Controller Design

The control objective is to design a robust control torqueτc such that the spacecraft can achieve the arbitrary attitude trajectory tracking demands in the presence of time-varying disturbances.For this reason, letQd=q0d,qdTT=q0d,q1d,q2d,q3dT∈ℝ4 be the desired quaternion and let ωd=ωxd,ωyd,ωzdT∈ℝ3 be the desired angular velocity. Define the tracking error vectors Qe=q0e,qeTT∈ℝ4 and ωe=ωxe,ωye,ωzeT∈ℝ3 as follows [36]:(70)Qe=Qd−1⊗Q=q0eqe=q0dq0+qdTqq0dq−q0qd−qd×q,ωe=ω−ωd=ωxeωyeωze=ωx−ωxdωy−ωydωz−ωzd,where ⊗ represents the quaternion multiplication. Based on (70), the quaternion-based error dynamics of (6), (7), and (25) is given by(71)q˙0e=−12qeTωe,(72)q˙e=12q0eI3+qe×ωe,(73)ω.e=Jeq−1−τc+d−ω×Jeqω+ΓJmΩ−ω.d.Select the sliding surface as(74)S=ωe+λqe,where S=s1,s2,s3T∈ℝ3; λ∈ℝ1 is a positive parameter to be designed.Suppose that the sliding motion is fulfilled in a finite timet=tf; it gives(75)S=S˙=0,∀t≥tf.From (71), (72), (74), and (75), the nonlinear reduced-order dynamics can be obtained:(76)q˙0e=−12qeTωe=−12qeT−λqe=12λqeTqe,∀t≥tf,(77)q˙e=12q0eI3+qe×−λqe=−12λq0eqe,∀t≥tf.The term “reduced-order” means that the system error dynamics described by (71)–(73 with order 7 “reduced” to subsystem (76) and (77) with order 4. In order to analyze the stability of the reduced-order dynamics, apply the identity of unit quaternion:(78)q0e2+qeTqe=1.Equation (76) can be decoupled as(79)q˙0e=12λ1−q0e2,∀t≥tf.Then, consider the following variable transformation:(80)q0et=1+1yt,which implies(81)q˙0et=−y˙ty2.Substituting (79) yields(82)−y˙ty2=12λ1−1+1yt2=12λ−2y−1y2,which implies(83)y˙−λy=12λ.Clearly, equation (83) is a linear ordinary differential equation. Its solution is given by(84)yt=Ceλt−12,where C is an integration constant. Apply the inverse mapping from (80); we get(85)q0et=1+1−1/2+Ceλt.Letq0e0 be the initial condition; we have(86)C=12+1q0e0−1.Hence, the analytic solution of (79) can be obtained:(87)q0et=1+2q0etf−1−q0etf+1+q0etf+1eλt−tf,∀t≥tf.Observing (87), it can be deduced that (i) the singularity occurs as an improper sliding gain λ is chosen such that −q0etf+1+q0etf+1eλt−tf=0,∀t≥tf. Thus, the sliding gain λ>0 is designed so that −q0etf+1+q0etf+1eλt−tf≠0,∀t≥tf. (ii) q0e⟶+1 as t⟶∞ from (87) and qe⟶0 as t⟶∞ from the identity of the unit quaternion (78). Observing (70), it can be found that q0e⟶1 and qe⟶0 imply q0⟶q0d and q⟶⟶q⟶d. That is, the nonlinear reduced-order dynamics (76) and (77) are asymptotically stable. It is different and more outstanding than [30, 35, 36]; the assumption that q0e≠0 is not made. In the following, how to enter the sliding mode in a finite time in the presence of external disturbances by means of the robust control law will be discussed.In order to introduceτc, taking the time derivative about (74) yields(88)S˙=ω.e+λq˙e=Jeq−1−τc+d−ω×Jeqω+ΓJmΩ+ξ,where(89)ξ=−ω.d+0.5λq0eI3+qe×ωe.Based on the super-twisting sliding mode algorithm [21], the following robust control law is designed:(90)τc0=ω×Jeqω+ΓJmΩ−Jeqξ,τcN=−JeqK1SS1/2+K2∫0tsignSτdτ,τc=−τc0+τcN,in which the gain matrices are denoted as(91)K1=diagk11k12k13,K2=diagk21k22k23.Substituting (90) into (88) yields the closed-loop sliding dynamics:(92)S˙=−K1SS1/2−K2∫0tsignSτdτ+D,where D=Jeq−1d=D1,D2,D3T∈ℝ3 and signS∈ℝ3 is a sign function defined as follows:(93)signsi=1,ifsi>0,−1,ifsi<0,signsi∈−1,1,ifsi=0,i=1,2,3..According to the stability criteria derived from Theorems1 and 2, the following gains are chosen:(94)k2i>δi,k1i>4k2i,where δi=supD˙i, i=1,2,3.


## 4. Spacecraft Reference Command Generation

### 4.1. Reference Quaternion Generation from TNB Frame

In this section, we are going to illustrate how to obtain the reference quaternion from the TNB frame. The procedures are as follows:(1)
Reference orbital trajectory: Consider the governing equations of the orbital motion [41]:(95)r¨=−μr3r.wherer=X,Y,ZT∈ℝ3 is the position vector and its magnitude is r=r; μ is the gravitational parameter. The reference orbital trajectory can be obtained by integrating (95) in a given initial condition.(2)
Construction of the TNB Frame: According to geometry kinematics, the TNB frame can be constructed by the trajectory information. Letv=r˙ be the velocity vector and let a=r¨ be the acceleration vector. The unit tangent vector can be computed by(96)e^t=vv.The unit normal vector can be obtained by(97)e^n=anan.wherean=a−at is the normal acceleration and at=a⋅e^te^t is the tangential acceleration. Based on the definition of TNB frame, the binormal vector is(98)e^b=e^t×e^n.(3)
Reference Quaternion: Based on the rotational property, the direction cosine matrix (DCM) can be composed bye^t, e^n, and e^b; that is,(99)Cv=e^t⋮e^n⋮e^b.Applying the quaternion kinematics, the reference quaternion is(100)q0d=121+C11v+C22v+C33v1/2,q1d=14q0C32v−C23v,q2d=14q0C13v−C31v,q3d=14q0C21v−C12v.


### 4.2. Reference Angular Velocity Generation

To generate a feasible command trajectory, the following property is proposed [35].Property 1.
The matrixEQ defined in (9) has the following properties:(101)ETQEQ=I3,ddtETQQ˙=ETQQ¨
By using Property1, from (8), the angular velocity and its time derivative can be expressed as(102)ω=2ETQQ˙,ω.=2ETQQ¨
Hence, the formula of the desired angular velocityωd associated with the desired quaternion is(103)ωd=2ETQdQ˙d,ω.d=2ETQdQ¨dwhere(104)EQd=−q1d−q2d−q3dq0d−q3dq2dq3dq0d−q1d−q2dq1dq0d.


## 4.1. Reference Quaternion Generation from TNB Frame

In this section, we are going to illustrate how to obtain the reference quaternion from the TNB frame. The procedures are as follows:(1)
Reference orbital trajectory: Consider the governing equations of the orbital motion [41]:(95)r¨=−μr3r.wherer=X,Y,ZT∈ℝ3 is the position vector and its magnitude is r=r; μ is the gravitational parameter. The reference orbital trajectory can be obtained by integrating (95) in a given initial condition.(2)
Construction of the TNB Frame: According to geometry kinematics, the TNB frame can be constructed by the trajectory information. Letv=r˙ be the velocity vector and let a=r¨ be the acceleration vector. The unit tangent vector can be computed by(96)e^t=vv.The unit normal vector can be obtained by(97)e^n=anan.wherean=a−at is the normal acceleration and at=a⋅e^te^t is the tangential acceleration. Based on the definition of TNB frame, the binormal vector is(98)e^b=e^t×e^n.(3)
Reference Quaternion: Based on the rotational property, the direction cosine matrix (DCM) can be composed bye^t, e^n, and e^b; that is,(99)Cv=e^t⋮e^n⋮e^b.Applying the quaternion kinematics, the reference quaternion is(100)q0d=121+C11v+C22v+C33v1/2,q1d=14q0C32v−C23v,q2d=14q0C13v−C31v,q3d=14q0C21v−C12v.


## 4.2. Reference Angular Velocity Generation

To generate a feasible command trajectory, the following property is proposed [35].Property 1.
The matrixEQ defined in (9) has the following properties:(101)ETQEQ=I3,ddtETQQ˙=ETQQ¨
By using Property1, from (8), the angular velocity and its time derivative can be expressed as(102)ω=2ETQQ˙,ω.=2ETQQ¨
Hence, the formula of the desired angular velocityωd associated with the desired quaternion is(103)ωd=2ETQdQ˙d,ω.d=2ETQdQ¨dwhere(104)EQd=−q1d−q2d−q3dq0d−q3dq2dq3dq0d−q1d−q2dq1dq0d.


## 5. Numerical Simulation

According to [42], the numerical data are considered and summarized as follows.(i)
Moment of inertia of the spacecraft and reaction wheels is as follows:(105)J=353−1.53282−1.5230,Jw=0.1260000.0630000.063N‐m.(ii)
The external disturbance is chosen as(106)dt=0.005sin0.8tcos0.5tcos0.3tTN−m.(iii)
Initial conditions are the following [41]: r0=8000,0,6000T (km); r˙0=0,7,0T (km/s); ω0=7,−8,−7T (degree/sec); Q0=0.999,0.017,−0.035,−0.026T; and gravitational parameter μ=3.987×105 (N-km2/kg).(iv)
The control gains are(107)K1=diag0.5470.5470.547,K2=diag0.050.050.05,λ=0.05.The sampling rate is 400 Hz. Total time span is 4.17 hours. The simulation results are shown in Figures4–7. For convenience, only the transition response is shown.Figure 4
Spacecraft 3D trajectory.Figure 5
Quaternion tracking error.Figure 6
Sliding surface.Figure 7
Evolution of reaction torque.With regard to Figure4, it is demonstrated that the spacecraft tracks the reference TNB attitude trajectory successfully. The attitude of the spacecraft with respect to the global frame XYZ is represented as xyz. The response of quaternion tracking errors is shown in Figure 5. It can verify the derivation of analytic solution (87) that the scalar component q0e⟶1 and the vector component qe⟶0 as t⟶∞ and they satisfy the constraint q0e+qeTqe=1. Figure 6 illustrates the response of the sliding variable. It can be found that the converging speed of the sliding variables is fast even though the spacecraft is under the environment with the time-varying external disturbances. The evolution of the reaction torques is shown in Figure 7. It can be seen that the control signal is smooth and there is no chattering phenomenon. We can deduce that the controller based on the super-twisting sliding algorithm can be realized in the spacecraft attitude control problem.


## 6. Conclusion

In this study, the attitude dynamics based on the redundant reaction wheels configuration of the spacecraft is derived. To achieve full degree of freedom attitude tracking control, the quaternion kinematics is introduced. For practical realization purpose, the nonsquare FDM is reformulated as a square, invertible matrix by means of solving a static optimization problem. The robust, continuous super-twisting sliding mode algorithm is adopted in the attitude controller design to guarantee robust performance in the presence of exogenous disturbances. Furthermore, the corresponding finite time stability based on the LMI is provided in the sense of Lyapunov. Asymptotic stability of the nonlinear reduced-order dynamics is proven by means of an analytic solution without imposing any limitation in quaternion. In a space mission with this scenario, the TNB command generation strategy is presented. Finally, the simulation results verify that the spacecraft can achieve the arbitrary attitude trajectory tracking demands in the presence of the time-varying external disturbances.



---

*Source: 6644033-2021-03-16.xml*