# Attitude-Tracking Control with Path Planning for Agile Satellite Using Double-Gimbal Control Moment Gyros

Authors: Peiling Cui; Fuyu Liu
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

In view of the issue of rapid attitude maneuver control of agile satellite, this paper presents an attitude-tracking control algorithm with path planning based on the improved genetic algorithm, adaptive backstepping control as well as sliding mode control. The satellite applies double gimbal control moment gyro as actuator and is subjected to the external disturbance and uncertain inertia properties. Firstly, considering the comprehensive mathematical model of the agile satellite and the double gimbal control moment gyro, an improved genetic algorithm is proposed to solve the attitude path-planning problem. The goal is to find an energy optimal path which satisfies certain maneuverability under the constraints of the input saturation, actuator saturation, slew rate limit and singularity measurement limit. Then, the adaptive backstepping control and sliding mode control are adopted in the design of the attitude-tracking controller to track accurately the desired path comprised of the satellite attitude quaternion and velocity. Finally, simulation results indicate the robustness and good tracking performance of the derived controller as well as its ability to avert the singularity of double gimbal control moment gyro.


---

## Body

## 1. Introduction

Rapid attitude maneuver control as one of the key technologies of agile satellite, has been extensively studied. Because of the superior properties such as large output torque and momentum storage, control moment gyro (CMG) is the ideal actuator for the rapid attitude maneuver control of agile satellite [1]. Compared with single gimbal CMG (SGCMG), double gimbal CMG (DGCMG) has many advantages, including the output control torque with two degrees of freedom, nearly spherical angular momentum envelop, good configuration efficiency, no striking singularity problem, low computation complexity of the steering law, and so on [2]. In the early 1970s, DGCMG has been successfully applied. For example, three orthogonally-mounted DGCMGs were used in the NASA’s Skylab [3]. And the International Station is controlled by four parallel-mounted DGCMGs [4].Many researches on satellite attitude control using CMG as actuator has been presented abundantly in the past decades [5–8]. The existing control systems are usually separated into the attitude control law and the CMG steering law. When the control torque is calculated by the attitude control law, the geometric singularity problem of the CMG is not taken into account. As a result, the singularity problem cannot be considered from a global perspective. An effective singularity avoidance method is to use a tracking controller to track the planned attitude trajectories accurately. And the planned attitude trajectories are based on the comprehensive mathematic model of agile satellite and CMG in the whole maneuvering process [9–11]. The attitude path planning problem of the agile satellite using SGCMGs is converted to an optimal control problem in [10]. And the Legendre pseudospectral method is used to solve the optimal control problem so that the comprehensive optimal attitude path is obtained in the aspect of the time and energy. Finally, a sliding mode controller is designed to track accurately the desired path. A feedback control with nominal inputs for agile satellites using SGCMGs is proposed in [11]. The feedforward control law using an energy optimal path planned by Fourier basis algorithm is used to solve the singularity problem of the SGCMGs. Moreover, a feedback control system is also designed in order to acquire the robustness against errors and disturbances. However, the attitude paths obtained in [10, 11] are local optimal [12]. The Legendre pseudospectral method and the Fourier basis algorithm are fundamentally gradient methods. The gradient method is a local method in that a locally optimal solution will generally be obtained. Whereas, the heuristic method is a global technique. In recent years, a path planning problem is introduced to new heuristic techniques such as fuzzy logic [13], ant colony optimization [14, 15], neural network [16], simulated annealing [17], genetic algorithm (GA) [18–21], and so on. Each method has its own strength over others in certain aspects. However, GA has attracted significant attention because of efficient use of large numbers of parallel processors, each of implementation for both continuous and discrete problems, no requirement for the continuity in response functions, and more robust solution generations to search for global or near global solutions.Sliding mode control is widely applied in spacecraft attitude control for its strong robustness against disturbances and uncertainties when the system states are sliding on the sliding surface [22, 23]. However, these design methods require the information on the bounds of the uncertainties/disturbances for the computation of the controller gains. Nevertheless, it is difficult to know the bounds of the uncertainties/disturbances in practical situation. Unlike these methods, nonlinear adaptive control methods do not require these bounds. They include the adaptive mechanism for estimating the uncertainties/disturbances [24]. Therefore, a variety of adaptive controllers are developed [25, 26]. Meanwhile, the adaptive robust spacecraft controllers based on the combination of sliding mode control and adaptive control are also designed [27, 28]. However, the methods mentioned above do not explicitly investigate the effects of constant disturbance toques on the attitude control. More recently, the backstepping control method is proposed to address the attitude control problem [24, 29, 30], due to its remarkable capability in designing cascaded systems, and it is well established that the spacecraft systems just satisfy the cascaded structure. Furthermore, the added integrator with backstepping control can improve the robustness of the system against the constant disturbance torques. For drawbacks of the simple linear backstepping controller, such as sluggish motion and excessive control input, the arctangent nonlinear tracking function is adopted in [29, 30]. Combined with the existing results in the literature, there are four advantages which make the adaptive backstepping sliding mode control attractive. Firstly, the adaptive backstepping sliding mode controller can be implemented easily and systematically by incorporating the backstepping design processes. Secondly, the convergence of the attitude tracking error for possible initial conditions can be guaranteed by selecting the Lyapunov function. Thirdly, the bounds of the uncertainties/disturbances do not require to be known because the uncertainties/disturbances can be estimated by the adaptation mechanism. Finally, the system robustness against modeling uncertainties and external disturbances can be improved because of the added integrator with backstepping control, so that the accuracy of steady-state control can be enhanced.This paper addresses the design of the closed-loop control system of the rapid attitude maneuver of agile satellite using DGCMGs. An improved GA is adopted to plan the energy optimal attitude trajectory based on the comprehensive model of agile satellite and DGCMGs. Moreover, the desired path planned by the improved GA must satisfy the variable constrained conditions, such as control input saturation, actuator saturation, slew rate limit, and singularity measurement limit. Then, an attitude tracking controller is proposed based on the sliding mode control method and adaptive backstepping control method using hyperbolic tangent function as a nonlinear tracking function. Finally, to realize the rapid attitude maneuver, the attitude tracking controller is applied to track accurately the desired path composed of the satellite attitude quaternion and velocity.


## 2. Dynamic Modeling of Rigid Agile Satellite Using DGCMGs

In this paper, the agile satellite has a DGCMG system that is an orthogonal configuration of two DGCMGs in Figure1. This section describes the dynamic model of agile satellite and the output torque equation of DGCMGs. The objective of this section is to construct a numerical calculation model of the satellite attitude maneuver using DGCMGs.Figure 1
Two DGCMGs orthogonal configuration.


### 2.1. Dynamic Modeling of Rigid Agile Satellite

The unit quaternion is adopted to describe the attitude of agile satellite to obtain the global representation without singularities. The unity quaternionQ is defined as Q=[q0q]T. Where q0 and q=[q1q2q3]T are the scalar and vector components of the unit quaternion, respectively. Let Qd=[qd0qd]T denote the desired quaternion. The error quaternion Qe=[qe0qe]T is described as
(2.1)qe0=q0qd0+qTqd,qe=qd0q-q0qd+q×qd.Letω∈R3 and ωd∈R3 denote the satellite angular velocity and the desired one. q× is the skew-symmetric matrix of the vector q. The skew-symmetric matrix of the vector a=[a1a2a3]T is defined by
(2.2)a×=[0-a3a2a30-a1-a2a10].Then the error of angular velocity is given as(2.3)ωe=ω-Cdbωd,
where Cdb=(qe02-qeTqe)I3×3+2qeqeT-2qe0qe×, I3×3 represents the identity matrix.The kinematic differential equations can be written as [31]
(2.4)q˙e0=-12qeTωe,q˙e=12Ξ(qe)ωe,
where
(2.5)Ξ(qe)=qe0I3×3+qe×.The control torque of DGCMG is generated through changing the direction of the angular momentum. Leth denote the total angular momentum of the DGCMGs, then the output torque of the DGCMGs can be represented as
(2.6)T=-h˙.According to the Euler equation, the error dynamic model can be expressed as [31]
(2.7)Jω˙e=Γ(J0,ωe)+u+Fd,
where
(2.8)Γ(J0,ωe)=-(ωe+Cdbωd)×[J0(ωe+Cdbωd)+h]+J0(ωe×Cdbωd-Cdbω˙d),Fd=-(ωe+Cdbωd)×[ΔJ(ωe+Cdbωd)]+ΔJ(ωe×Cdbωd-Cdbω˙d)+Td,
where J=J0+ΔJ is the satellite inertia matrix, and J0 and ΔJ are the nominal and uncertain components of the inertia matrix. u=T∈R3 is the command torque, Td∈R3 is the external disturbance torque, including the gravity gradient torque, aerodynamic torque, and torque of solar radiation. It is assumed that Td is an unknown bounded function. Fd that results from ΔJ and Td is the uncertain part in the dynamic equation, and it is also an unknown bounded function.


### 2.2. Output Torque Equation of DGCMGs

In order to realize the three-axis attitude control, an orthogonal configuration of two DGCMGs is adopted. The configuration diagram is given in Figure1.The vector set{Xb,Yb,Zb} is the unit vector of the satellite body coordinate frame. Xb, Yb, and Zb are orthogonal to each other. xi, yi, and zi are parallel to the outer gimbal axis, spin axis and inner gimbal axis of the ith DGCMG, i=1,2. αi and βi are the inner gimbal angle and outer gimbal angle of the ith DGCMG. Let h0 denote the angular momentum magnitude of the individual DGCMG. Then the angular momentum of the DGCMGs is a function depending on the gimbal angles δ=[δ1δ2]T. It is expressed as
(2.9)h=h0[cosα1cosβ1+cosα2cosβ2-cosα1sinβ1-sinα2-sinα1+cosα2sinβ2],
where δi=[αiβi]T, i=1,2. According to the principle of angular momentum exchange, the time derivative of h is equal and opposite to the output torque T from the DGCMGs to the satellite body, and their relationship is
(2.10)T=-h˙=-h0Cδ˙,
where
(2.11)C=[-sinα1cosβ1-cosα1sinβ1-sinα2cosβ2-cosα2sinβ2sinα1sinβ1-cosα1cosβ1-cosα20-cosα10-sinα2sinβ2cosα2cosβ2].An obstacle when using a CMG system in practice is the existence of singular gimbal angle states for which the CMGs cannot generate torque along arbitrary directions. The singularity measurement is given as [4]
(2.12)d=det(CCT).This determinant is an indicator of how close the CMG cluster is to a singular point. Whered is a positive scalar. When d=0, DGCMGs are strapped into the singularity state. The distance from the singularity state is farther when d is larger.


## 2.1. Dynamic Modeling of Rigid Agile Satellite

The unit quaternion is adopted to describe the attitude of agile satellite to obtain the global representation without singularities. The unity quaternionQ is defined as Q=[q0q]T. Where q0 and q=[q1q2q3]T are the scalar and vector components of the unit quaternion, respectively. Let Qd=[qd0qd]T denote the desired quaternion. The error quaternion Qe=[qe0qe]T is described as
(2.1)qe0=q0qd0+qTqd,qe=qd0q-q0qd+q×qd.Letω∈R3 and ωd∈R3 denote the satellite angular velocity and the desired one. q× is the skew-symmetric matrix of the vector q. The skew-symmetric matrix of the vector a=[a1a2a3]T is defined by
(2.2)a×=[0-a3a2a30-a1-a2a10].Then the error of angular velocity is given as(2.3)ωe=ω-Cdbωd,
where Cdb=(qe02-qeTqe)I3×3+2qeqeT-2qe0qe×, I3×3 represents the identity matrix.The kinematic differential equations can be written as [31]
(2.4)q˙e0=-12qeTωe,q˙e=12Ξ(qe)ωe,
where
(2.5)Ξ(qe)=qe0I3×3+qe×.The control torque of DGCMG is generated through changing the direction of the angular momentum. Leth denote the total angular momentum of the DGCMGs, then the output torque of the DGCMGs can be represented as
(2.6)T=-h˙.According to the Euler equation, the error dynamic model can be expressed as [31]
(2.7)Jω˙e=Γ(J0,ωe)+u+Fd,
where
(2.8)Γ(J0,ωe)=-(ωe+Cdbωd)×[J0(ωe+Cdbωd)+h]+J0(ωe×Cdbωd-Cdbω˙d),Fd=-(ωe+Cdbωd)×[ΔJ(ωe+Cdbωd)]+ΔJ(ωe×Cdbωd-Cdbω˙d)+Td,
where J=J0+ΔJ is the satellite inertia matrix, and J0 and ΔJ are the nominal and uncertain components of the inertia matrix. u=T∈R3 is the command torque, Td∈R3 is the external disturbance torque, including the gravity gradient torque, aerodynamic torque, and torque of solar radiation. It is assumed that Td is an unknown bounded function. Fd that results from ΔJ and Td is the uncertain part in the dynamic equation, and it is also an unknown bounded function.


## 2.2. Output Torque Equation of DGCMGs

In order to realize the three-axis attitude control, an orthogonal configuration of two DGCMGs is adopted. The configuration diagram is given in Figure1.The vector set{Xb,Yb,Zb} is the unit vector of the satellite body coordinate frame. Xb, Yb, and Zb are orthogonal to each other. xi, yi, and zi are parallel to the outer gimbal axis, spin axis and inner gimbal axis of the ith DGCMG, i=1,2. αi and βi are the inner gimbal angle and outer gimbal angle of the ith DGCMG. Let h0 denote the angular momentum magnitude of the individual DGCMG. Then the angular momentum of the DGCMGs is a function depending on the gimbal angles δ=[δ1δ2]T. It is expressed as
(2.9)h=h0[cosα1cosβ1+cosα2cosβ2-cosα1sinβ1-sinα2-sinα1+cosα2sinβ2],
where δi=[αiβi]T, i=1,2. According to the principle of angular momentum exchange, the time derivative of h is equal and opposite to the output torque T from the DGCMGs to the satellite body, and their relationship is
(2.10)T=-h˙=-h0Cδ˙,
where
(2.11)C=[-sinα1cosβ1-cosα1sinβ1-sinα2cosβ2-cosα2sinβ2sinα1sinβ1-cosα1cosβ1-cosα20-cosα10-sinα2sinβ2cosα2cosβ2].An obstacle when using a CMG system in practice is the existence of singular gimbal angle states for which the CMGs cannot generate torque along arbitrary directions. The singularity measurement is given as [4]
(2.12)d=det(CCT).This determinant is an indicator of how close the CMG cluster is to a singular point. Whered is a positive scalar. When d=0, DGCMGs are strapped into the singularity state. The distance from the singularity state is farther when d is larger.


## 3. Attitude Tracking Control with Path-Planning

The whole attitude control system consists of the attitude path planning and the attitude tracking control, which applies DGCMGs as actuator. In the attitude path planning, an improved GA algorithm is adopted to find the energy optimal path without separating the whole control system into the attitude control law and the CMG steering law. Moreover, the desired path planned by GA satisfies certain maneuver ability and various physical constrains, including input saturation, actuator saturation, angular velocity limit and singularity measurement limit. In the attitude tracking control, to guarantee the system robustness to the external disturbance and uncertain inertia properties, an adaptive backstepping sliding mode controller is designed to track accurately the desired path comprised of the attitude quaternion and velocity instead of gimbal angular velocity of the DGCMGs. Using the attitude quaternion and velocity as reference trajectory will lead to the few deviation of gimbal angular velocity of the DGCMGs from the desired path in the practical maneuvering process. Consequently, the singularity of DGCMGs cannot be avoided completely. However, the robustness of the system will be enhanced. The whole control system block diagram is shown in Figure2.Figure 2
Block diagram of the proposed control system.


### 3.1. Reference Maneuver Path Planning Using GA

This section deals with the trajectory planning problem for agile satellite using DGCMGs as actuator. An improved GA using a floating representation is proposed to search for the energy optimal path. To overcome the premature convergence or get into local optimal problem of the classic GA, considering the characteristics of rapid attitude maneuver of the agile satellite, a heuristic method is utilized to produce the initial population. Then, in order to guarantee the DGCMGs away from the singularity state, a fitness function is defined based on the energy consumed by DGCMGs. Moreover, for purpose of making the command torque and the angular velocity of satellite smoothing, the smooth mutation operator for the angular velocity and the angular acceleration is presented.(1) Individual Representation
Compared with the binary representation, the floating representation can be implemented easily, and computer memory is saved. Moreover, it is more suitable for solving multidimensional, high-precision numerical problem [18]. Therefore, the floating representation is used to encode attitude quaternion path. The chromosome structure is denoted as (qe,1,qe,2,…,qe,n). Where qe,1 and qe,n are the vector components of the initial attitude and target attitude, respectively. qe,i(i=1,2,…n-1) is the vector components of the attitude quaternion of the intermediate nodes. The path length n is calculated according to the simulation time T and the simulation step Ts, n=T/Ts.(2) Population Initialization
In order to improve the searching efficiency, a heuristic method is adopted to produce the initial population. In this problem, the boundary conditions of the state variables are described as(3.1)qe(0)=qe0,qe(T)=qed,ω(0)=ω0,ω(T)=ωd,ω˙(T)=ω˙d.To satisfy the boundary conditions,qe,2, qe,n-2, and qe,n-1 are calculated based on the comprehensive mathematic model of the agile satellite using DGCMGs as actuator. In addition, qe,i(i=3,…n-3) is produced by the following rules:(1)
To make the energy minimal,qe→0 is monotone increasing or decreasing. So part of the initial population is generated as follows:
(3.2)qe,i=(qe,i-1-qe,n-2)x+qe,n-2,x∈(0,1),i∈(3,4,…,n-3).(2)
If the singular point is not taken into account, the connection betweenqe,2 and qe,n-2 is the shortest and the most smooth. Based on this idea, part of initial population is produced.(3)
When the agile satellite realizes the fastest maneuver, the angular velocity of the agile satellite usually includes acceleration and deceleration that are similar to theS function derivative curve. Therefore, part of the initial population is initialized as follows:
(3.3)qe,i=sign(qe,2)|qe,2|1+expn[-(ciTs+(a-2)Ts-m)]+qe,2,i∈(3,…,n-3),
where m∈(0,T/2), n∈(0,30). a and c are given as
(3.4)f2(t)=sign(qe,2)|qe,2|1+expn[-(t-m)]+qe,2,Δf2(i)=f2((i+1)Ts)-f2(iTs),a=min(find(Δf2>1e-3)),b=max(find(Δf2>1e-3)),c=b-a+1n-4
where t∈(0,Ta), Ta is an appropriate constant greater than T.(3) Fitness Function
A fitness function is designed based on the energy consumed by DGCMGs for ensuring the DGCMGs away from the singularity state. Usually, the energy consumption is estimated by the norm of control inputs. However, the energy consumption of DGCMGs cannot be estimated precisely because the satellite angular velocities cannot be ignored with respect to the DGCMGs gimbal angular velocities in the agile satellite. Therefore, the energy consumption of the DGCMGs required to run a given pathIj is set up as follows:
(3.5)E(Ij)=∑k=1nWk′Wk,Wk=Tk-ωk×hk,
where Wk is an output torque vector of the DGCMGs on the kth point in the path Ij. This evaluation term for the energy consumption is derived from an assumption of the proportional property between the output torque and the electrical energy consumption of the DGCMG gimbal [11].Moreover, the trajectory planned by the improved GA must satisfy various physical constrained conditions, including the input saturation, gimbal angular velocity of DGCMGs limit, angular velocity of agile satellite limit, and singularity measurement limit. Then the evaluation function is adjusted to the following constrained conditions:(3.6)Ω1={u∣|ui|<um,i=1,2,3},Ω2={δ˙∣|δ˙i|<δ˙m,i=1,2,3,4},Ω3={ω∣|ωi|<ωm,i=1,2,3},d>dmin,
where um, δ˙m, ωm and dmin denote the maximal command torque, the maximal gimbal angular velocities of the DGCMGs, the maximal angular velocity of satellite, and the minimal singularity measurement, respectively. To guarantee the system robustness, um, δ˙m, and ωm are less than the actual maximal command umax, actual maximal gimbal angular velocities of DGCMGs δ˙max, and the actual maximal angular velocity of satellite ωmax, respectively.To satisfy these constrained conditions with keeping the energy consumption of the DGCMGs low, the evaluation function of this optimization is modified in the following form:(3.7)F(Ij)=∑k=1nWk′Wk+num×inf1,Wk=Tk-ωk×hk,
where num is the number of the unfeasible points in the path Ij, and inf1 is an appropriate positive constant. Equation (3.7) represents that those paths with unfeasible nodes are heavily penalized with an extra energy cost.(4) Next Generation Process
For purpose of making the command torque and the angular velocity of satellite smoothing, a smooth mutation operator for the angular velocity and angular acceleration is presented.Selection. The N/2 individuals with the lowest energy cost (strongest) are chosen from the current population.Crossover. New individuals (called offspring) are generated from the selected path by using a crossover operator. Therefore, a total of N/4 pairs are formed. The first two offspring of each pair are equal to their parents. The two other offspring are produced as the recombination of their parents. The recombination process between the two selected paths is defined as follows [18].Supposing the selected two chromosomes are(3.8)Ii={qe,1,iqe,2,i…qe,(n-1),iqe,n,i},Ij={qe,1,jqe,2,j…qe,(n-1),jqe,n,j}.Thus supposing the split position is ati, the offspring is produced as follows:
(3.9)Ii′={(1-α)qe,1,i+αqe,1,j…(1-α)qe,n,i+αqe,n,j},Ij′={(1-α)qe,1,j+αqe,1,i…(1-α)qe,n,j+αqe,n,i},
where α∈(0,1). The greatest possible information exchange is obtained by choosing α that is very important for the crossover operator. In this paper, α is randomly selected. However, for different chromosomes, α is different.Mutation. To satisfy the boundary conditions described as (3.1), the mutation operation is performed on qe,i(i=3,…n-3). When the path is feasible, a small percentage of feasible points is mutated at random and ensure that the path is still feasible. On the other hand, when the path is unfeasible, a large percentage of unfeasible points is mutated at random.Moreover, due to the smooth requirement of the angular velocity and the command torque, the smooth mutation operator for the angular velocity and the angular acceleration is executed. If one point that is more or less than its adjacent points is found, this point is mutated between the two adjacent points.Finally, to avoid that the strongest individual is destroyed, the excellent individual protection method is adopted. The best individual of each generation is directly inherited into the next generation, and not be mutated. This method will improve the convergence speed of the algorithm.


### 3.2. Attitude Tracking Controller Design

In the attitude path planning, various physical constrained conditions are considered and have a certain margin for the external disturbances and uncertainties. Therefore, these physical constrained conditions including the input limitation is not required to be considered in the attitude tracking controller design. In order to track accurately the desired path, and to overcome the external disturbance and the uncertain inertial properties, the attitude tracking controller is derived based on the sliding mode control method and the backstepping control method using hyperbolic tangent function as a nonlinear tracking function. Because it is difficult to know the uncertain inertial and the external disturbances in practical situation, it is necessary to assume that the uncertainties that result from the external disturbance and inertial properties is an unknown bounded function. Then, the adaptive control method is utilized to estimate the uncertainties that result from the external disturbance and inertial properties.(1) Adaptive Backstepping Sliding Mode Controller Design Method
The attitude control system described by (2.4) and (2.7) has a cascaded structure. The error of the angular velocity vector ωe can be thought of as a pseudoinput vector, while the command torque u is a real control input vector that affects ωe. Thus, the backstepping control method can be easily employed in the satellite attitude control system.The first step in the backstepping control method is to find the control law for the subsystem described by (2.4). The state variable x1 is defined by
(3.10)x1=qe.Let us consider the following candidate Lyapunov function(3.11)V1=12x1Tx1+12(1-qe0)2.Substituting (2.4) into the time derivative of (3.11) yields
(3.12)V˙1=12x1Tωe.In accordance with the discussion performed in [29], the simple linear tracking law is not effective for the satellite maneuver for the reason of sluggish motion, trivial nonlinear term cancellation, and excessive control input  at the initial stage of the maneuver. To overcome these defects, an effective nonlinear tracking function ϕ(x1) is proposed. When it is adopted, the nonlinear tracking law is obtained:
(3.13)ωe=-kϕ(x1).Supporting the nonlinear functionϕ(x1) satisfies
(3.14)V˙1=-12kx1Tϕ(x1)≤0,
where k>0. Only when x1=0, V˙1=0. Then x1 is asymptotic stability.After designing the tracking law, the real control inputu should be determined such that (3.13) is achieved. The state variable x2 is selected as
(3.15)x2=ωe+kϕ(x1).
By choosing the following sliding surface
(3.16)s=cx1+x2,
where c>0, then
(3.17)Js˙=-12kJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)+12J0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2+Γ(J0,ωe)+u+F,
where ∂ϕ(x1)/∂x1 is the Jacobian matrix of ϕ(x1), F is the total uncertain, and is expressed as
(3.18)F=Fd-12kΔJΞ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)+12ΔJΞ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2.BecauseΔJ and Td are difficult to be known in practical situation, the uncertainties F that result from ΔJ and Td is an unknown bounded function. Consider the following candidate Lyapunov function for the overall system
(3.19)V2=V1+12sTJs+12γF~TF~,
where γ is a positive constant. F^ is the estimation of F, F~=F-F^ is the estimation error of F. Because the uncertain inertial momentum and external disturbance are changing gradually, F~˙=-F^˙.The derivative of (3.19) is obtained
(3.20)V˙2=-12ksTJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)+12sTJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2+12x1Tx2+sT(Γ(J0,ωe)+u+F^)-12kx1Tϕ(x1)-1γF~T(F^˙-γs).Letting(3.21)u=12kJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)-Γ(J0,ωe)-F^-12J0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2-g(s+βsign(s)).For the reason thatF is an unknown bounded function, an adaptive law is utilized to estimate it. And the adaptive law is selected as follows:
(3.22)F^˙=γs,
Where g, β and γ are positive constants. Substituting (3.21) and (3.22) into (3.20) gives
(3.23)V˙2=12x1Tx2-12kx1Tϕ(x1)-gsTs-gβsTsign(s).Letting(3.24)R=[12k+gc2gc-14gc-14g]Owing to(3.25)xTRx=-12x1Tx2+12kx1Tx1+gsTs,
Where x=[x1x2]T. Then (3.23) can be rewritten as
(3.26)V˙2=-xTRx-12kx1T(ϕ(x1)-x1)-gβsTsign(s).Duing to(3.27)|R|=(12k+gc2)g-(gc-14)2=12g(k+c)-116.By choosing the values ofg, c, and k, it can be guaranteed that R is a positive definite matrix. Simultaneously, x1T(ϕ(x1)-x1)≥0 can be guaranteed by choosing the nonlinear tracking function ϕ(x1). Therefore, V˙2≤0 can be ensured by selecting the values of  g, c, k, β, and the nonlinear tracking function ϕ(x1).(2) Nonlinear Tracking Function Selection
According to the design process of the adaptive backstepping sliding mode controller, to guaranteeV˙1≤0 and V˙2≤0, the nonlinear tracking function ϕ(x1) must satisfy
(3.28)x1Tϕ(x1)≥0,x1T(ϕ(x1)-x1)≥0.According to [29], the following condition can contribute to reduce the peak control torque
(3.29)ϕ(x1)≈1,ifqei=[1-μ,1],ϕ(x1)≈-1,ifqei=[-1,-1+μ],
where 0<μ<1. The above condition prescribes that the nonlinear tracking function value does not produce a large value when |qei|→1. Meanwhile, if the following condition is imposed, the convergence speed can be increased when |qei|→0. 								(3.30)ϕ(x1)>κx1,ifqei→+0,ϕ(x1)<-κx1,ifqei→-0,
where κ>1. Equation (3.30) requires that ϕ(x1) produce higher absolute value than the linear tracking function when |qei| is small. This condition contributes to the removal of sluggish motion. The following hyperbolic tangent function would be one of the candidate functions that satisfy (3.28)–(3.30)
(3.31)ϕ(x1)=tanh(αx1),
where α≥4. In this paper, (3.31) is chosen as the nonlinear tracking function because of its simplicity and smoothness. Moreover, it has fewer parameters than the arctangent function. Differentiating (3.31) with respect to time gives
(3.32)∂ϕ(x1)∂x1=diag{αsech2(qe1),αsech2(qe2),αsech2(qe3)}.To avoid the chattering of the control input, the sign function in (3.21) can be replaced by
(3.33)sat(si,ξ)={1si>ξ,siξ|si|≤ξ,-1si<-ξ,
where ξ is a small positive constant. Equation (3.21) can be rewritten as
(3.34)u=12kJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)-Γ(J0,ωe)-F^-12J0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2-g(s+βsat(s)).From above, the adaptive backstepping sliding mode controller includes (3.22), (3.31), (3.32) and (3.34).


## 3.1. Reference Maneuver Path Planning Using GA

This section deals with the trajectory planning problem for agile satellite using DGCMGs as actuator. An improved GA using a floating representation is proposed to search for the energy optimal path. To overcome the premature convergence or get into local optimal problem of the classic GA, considering the characteristics of rapid attitude maneuver of the agile satellite, a heuristic method is utilized to produce the initial population. Then, in order to guarantee the DGCMGs away from the singularity state, a fitness function is defined based on the energy consumed by DGCMGs. Moreover, for purpose of making the command torque and the angular velocity of satellite smoothing, the smooth mutation operator for the angular velocity and the angular acceleration is presented.(1) Individual Representation
Compared with the binary representation, the floating representation can be implemented easily, and computer memory is saved. Moreover, it is more suitable for solving multidimensional, high-precision numerical problem [18]. Therefore, the floating representation is used to encode attitude quaternion path. The chromosome structure is denoted as (qe,1,qe,2,…,qe,n). Where qe,1 and qe,n are the vector components of the initial attitude and target attitude, respectively. qe,i(i=1,2,…n-1) is the vector components of the attitude quaternion of the intermediate nodes. The path length n is calculated according to the simulation time T and the simulation step Ts, n=T/Ts.(2) Population Initialization
In order to improve the searching efficiency, a heuristic method is adopted to produce the initial population. In this problem, the boundary conditions of the state variables are described as(3.1)qe(0)=qe0,qe(T)=qed,ω(0)=ω0,ω(T)=ωd,ω˙(T)=ω˙d.To satisfy the boundary conditions,qe,2, qe,n-2, and qe,n-1 are calculated based on the comprehensive mathematic model of the agile satellite using DGCMGs as actuator. In addition, qe,i(i=3,…n-3) is produced by the following rules:(1)
To make the energy minimal,qe→0 is monotone increasing or decreasing. So part of the initial population is generated as follows:
(3.2)qe,i=(qe,i-1-qe,n-2)x+qe,n-2,x∈(0,1),i∈(3,4,…,n-3).(2)
If the singular point is not taken into account, the connection betweenqe,2 and qe,n-2 is the shortest and the most smooth. Based on this idea, part of initial population is produced.(3)
When the agile satellite realizes the fastest maneuver, the angular velocity of the agile satellite usually includes acceleration and deceleration that are similar to theS function derivative curve. Therefore, part of the initial population is initialized as follows:
(3.3)qe,i=sign(qe,2)|qe,2|1+expn[-(ciTs+(a-2)Ts-m)]+qe,2,i∈(3,…,n-3),
where m∈(0,T/2), n∈(0,30). a and c are given as
(3.4)f2(t)=sign(qe,2)|qe,2|1+expn[-(t-m)]+qe,2,Δf2(i)=f2((i+1)Ts)-f2(iTs),a=min(find(Δf2>1e-3)),b=max(find(Δf2>1e-3)),c=b-a+1n-4
where t∈(0,Ta), Ta is an appropriate constant greater than T.(3) Fitness Function
A fitness function is designed based on the energy consumed by DGCMGs for ensuring the DGCMGs away from the singularity state. Usually, the energy consumption is estimated by the norm of control inputs. However, the energy consumption of DGCMGs cannot be estimated precisely because the satellite angular velocities cannot be ignored with respect to the DGCMGs gimbal angular velocities in the agile satellite. Therefore, the energy consumption of the DGCMGs required to run a given pathIj is set up as follows:
(3.5)E(Ij)=∑k=1nWk′Wk,Wk=Tk-ωk×hk,
where Wk is an output torque vector of the DGCMGs on the kth point in the path Ij. This evaluation term for the energy consumption is derived from an assumption of the proportional property between the output torque and the electrical energy consumption of the DGCMG gimbal [11].Moreover, the trajectory planned by the improved GA must satisfy various physical constrained conditions, including the input saturation, gimbal angular velocity of DGCMGs limit, angular velocity of agile satellite limit, and singularity measurement limit. Then the evaluation function is adjusted to the following constrained conditions:(3.6)Ω1={u∣|ui|<um,i=1,2,3},Ω2={δ˙∣|δ˙i|<δ˙m,i=1,2,3,4},Ω3={ω∣|ωi|<ωm,i=1,2,3},d>dmin,
where um, δ˙m, ωm and dmin denote the maximal command torque, the maximal gimbal angular velocities of the DGCMGs, the maximal angular velocity of satellite, and the minimal singularity measurement, respectively. To guarantee the system robustness, um, δ˙m, and ωm are less than the actual maximal command umax, actual maximal gimbal angular velocities of DGCMGs δ˙max, and the actual maximal angular velocity of satellite ωmax, respectively.To satisfy these constrained conditions with keeping the energy consumption of the DGCMGs low, the evaluation function of this optimization is modified in the following form:(3.7)F(Ij)=∑k=1nWk′Wk+num×inf1,Wk=Tk-ωk×hk,
where num is the number of the unfeasible points in the path Ij, and inf1 is an appropriate positive constant. Equation (3.7) represents that those paths with unfeasible nodes are heavily penalized with an extra energy cost.(4) Next Generation Process
For purpose of making the command torque and the angular velocity of satellite smoothing, a smooth mutation operator for the angular velocity and angular acceleration is presented.Selection. The N/2 individuals with the lowest energy cost (strongest) are chosen from the current population.Crossover. New individuals (called offspring) are generated from the selected path by using a crossover operator. Therefore, a total of N/4 pairs are formed. The first two offspring of each pair are equal to their parents. The two other offspring are produced as the recombination of their parents. The recombination process between the two selected paths is defined as follows [18].Supposing the selected two chromosomes are(3.8)Ii={qe,1,iqe,2,i…qe,(n-1),iqe,n,i},Ij={qe,1,jqe,2,j…qe,(n-1),jqe,n,j}.Thus supposing the split position is ati, the offspring is produced as follows:
(3.9)Ii′={(1-α)qe,1,i+αqe,1,j…(1-α)qe,n,i+αqe,n,j},Ij′={(1-α)qe,1,j+αqe,1,i…(1-α)qe,n,j+αqe,n,i},
where α∈(0,1). The greatest possible information exchange is obtained by choosing α that is very important for the crossover operator. In this paper, α is randomly selected. However, for different chromosomes, α is different.Mutation. To satisfy the boundary conditions described as (3.1), the mutation operation is performed on qe,i(i=3,…n-3). When the path is feasible, a small percentage of feasible points is mutated at random and ensure that the path is still feasible. On the other hand, when the path is unfeasible, a large percentage of unfeasible points is mutated at random.Moreover, due to the smooth requirement of the angular velocity and the command torque, the smooth mutation operator for the angular velocity and the angular acceleration is executed. If one point that is more or less than its adjacent points is found, this point is mutated between the two adjacent points.Finally, to avoid that the strongest individual is destroyed, the excellent individual protection method is adopted. The best individual of each generation is directly inherited into the next generation, and not be mutated. This method will improve the convergence speed of the algorithm.


## 3.2. Attitude Tracking Controller Design

In the attitude path planning, various physical constrained conditions are considered and have a certain margin for the external disturbances and uncertainties. Therefore, these physical constrained conditions including the input limitation is not required to be considered in the attitude tracking controller design. In order to track accurately the desired path, and to overcome the external disturbance and the uncertain inertial properties, the attitude tracking controller is derived based on the sliding mode control method and the backstepping control method using hyperbolic tangent function as a nonlinear tracking function. Because it is difficult to know the uncertain inertial and the external disturbances in practical situation, it is necessary to assume that the uncertainties that result from the external disturbance and inertial properties is an unknown bounded function. Then, the adaptive control method is utilized to estimate the uncertainties that result from the external disturbance and inertial properties.(1) Adaptive Backstepping Sliding Mode Controller Design Method
The attitude control system described by (2.4) and (2.7) has a cascaded structure. The error of the angular velocity vector ωe can be thought of as a pseudoinput vector, while the command torque u is a real control input vector that affects ωe. Thus, the backstepping control method can be easily employed in the satellite attitude control system.The first step in the backstepping control method is to find the control law for the subsystem described by (2.4). The state variable x1 is defined by
(3.10)x1=qe.Let us consider the following candidate Lyapunov function(3.11)V1=12x1Tx1+12(1-qe0)2.Substituting (2.4) into the time derivative of (3.11) yields
(3.12)V˙1=12x1Tωe.In accordance with the discussion performed in [29], the simple linear tracking law is not effective for the satellite maneuver for the reason of sluggish motion, trivial nonlinear term cancellation, and excessive control input  at the initial stage of the maneuver. To overcome these defects, an effective nonlinear tracking function ϕ(x1) is proposed. When it is adopted, the nonlinear tracking law is obtained:
(3.13)ωe=-kϕ(x1).Supporting the nonlinear functionϕ(x1) satisfies
(3.14)V˙1=-12kx1Tϕ(x1)≤0,
where k>0. Only when x1=0, V˙1=0. Then x1 is asymptotic stability.After designing the tracking law, the real control inputu should be determined such that (3.13) is achieved. The state variable x2 is selected as
(3.15)x2=ωe+kϕ(x1).
By choosing the following sliding surface
(3.16)s=cx1+x2,
where c>0, then
(3.17)Js˙=-12kJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)+12J0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2+Γ(J0,ωe)+u+F,
where ∂ϕ(x1)/∂x1 is the Jacobian matrix of ϕ(x1), F is the total uncertain, and is expressed as
(3.18)F=Fd-12kΔJΞ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)+12ΔJΞ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2.BecauseΔJ and Td are difficult to be known in practical situation, the uncertainties F that result from ΔJ and Td is an unknown bounded function. Consider the following candidate Lyapunov function for the overall system
(3.19)V2=V1+12sTJs+12γF~TF~,
where γ is a positive constant. F^ is the estimation of F, F~=F-F^ is the estimation error of F. Because the uncertain inertial momentum and external disturbance are changing gradually, F~˙=-F^˙.The derivative of (3.19) is obtained
(3.20)V˙2=-12ksTJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)+12sTJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2+12x1Tx2+sT(Γ(J0,ωe)+u+F^)-12kx1Tϕ(x1)-1γF~T(F^˙-γs).Letting(3.21)u=12kJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)-Γ(J0,ωe)-F^-12J0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2-g(s+βsign(s)).For the reason thatF is an unknown bounded function, an adaptive law is utilized to estimate it. And the adaptive law is selected as follows:
(3.22)F^˙=γs,
Where g, β and γ are positive constants. Substituting (3.21) and (3.22) into (3.20) gives
(3.23)V˙2=12x1Tx2-12kx1Tϕ(x1)-gsTs-gβsTsign(s).Letting(3.24)R=[12k+gc2gc-14gc-14g]Owing to(3.25)xTRx=-12x1Tx2+12kx1Tx1+gsTs,
Where x=[x1x2]T. Then (3.23) can be rewritten as
(3.26)V˙2=-xTRx-12kx1T(ϕ(x1)-x1)-gβsTsign(s).Duing to(3.27)|R|=(12k+gc2)g-(gc-14)2=12g(k+c)-116.By choosing the values ofg, c, and k, it can be guaranteed that R is a positive definite matrix. Simultaneously, x1T(ϕ(x1)-x1)≥0 can be guaranteed by choosing the nonlinear tracking function ϕ(x1). Therefore, V˙2≤0 can be ensured by selecting the values of  g, c, k, β, and the nonlinear tracking function ϕ(x1).(2) Nonlinear Tracking Function Selection
According to the design process of the adaptive backstepping sliding mode controller, to guaranteeV˙1≤0 and V˙2≤0, the nonlinear tracking function ϕ(x1) must satisfy
(3.28)x1Tϕ(x1)≥0,x1T(ϕ(x1)-x1)≥0.According to [29], the following condition can contribute to reduce the peak control torque
(3.29)ϕ(x1)≈1,ifqei=[1-μ,1],ϕ(x1)≈-1,ifqei=[-1,-1+μ],
where 0<μ<1. The above condition prescribes that the nonlinear tracking function value does not produce a large value when |qei|→1. Meanwhile, if the following condition is imposed, the convergence speed can be increased when |qei|→0. 								(3.30)ϕ(x1)>κx1,ifqei→+0,ϕ(x1)<-κx1,ifqei→-0,
where κ>1. Equation (3.30) requires that ϕ(x1) produce higher absolute value than the linear tracking function when |qei| is small. This condition contributes to the removal of sluggish motion. The following hyperbolic tangent function would be one of the candidate functions that satisfy (3.28)–(3.30)
(3.31)ϕ(x1)=tanh(αx1),
where α≥4. In this paper, (3.31) is chosen as the nonlinear tracking function because of its simplicity and smoothness. Moreover, it has fewer parameters than the arctangent function. Differentiating (3.31) with respect to time gives
(3.32)∂ϕ(x1)∂x1=diag{αsech2(qe1),αsech2(qe2),αsech2(qe3)}.To avoid the chattering of the control input, the sign function in (3.21) can be replaced by
(3.33)sat(si,ξ)={1si>ξ,siξ|si|≤ξ,-1si<-ξ,
where ξ is a small positive constant. Equation (3.21) can be rewritten as
(3.34)u=12kJ0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)ϕ(x1)-Γ(J0,ωe)-F^-12J0Ξ(qe)(cI3×3+k∂ϕ(x1)∂x1)x2-g(s+βsat(s)).From above, the adaptive backstepping sliding mode controller includes (3.22), (3.31), (3.32) and (3.34).


## 4. Numerical Simulation

This section presents the simulation results to illustrate the effectiveness of the attitude tracking control with path-planning.In the simulation, the agile satellite is assumed to be a rigid body, and two DGCMGs are installed as illustrated in Figure1. Consider the spacecraft model (2.4) and (2.7) with the nominal inertia matrix
(4.1)J=[120.120.120.13130.130.150.1515]kg·m2.The uncertainties of the inertia matrix is(4.2)ΔJ=[10.120.120.131.50.130.150.152]kg·m2.The other parameters for satellite and DGCMGs are given in Table1. The external disturbance is
(4.3)Td=[3cosω0t+21.5sinω0t+3cosω0t+1-3cosω0t+1]×10-3N·m,
where ω0=0.0011 rad/s. The parameters for attitude tracking controller are chosen as α=8, k=1.5, c=1.5, g=2.6, β=0.1, and γ=0.1.Table 1
Simulation parameters.
Simulation parameters
Value
angular momentum magnitude ofa DGCMG h0
6 N·m·s
Initial attitude angle[φ0θ0ψ0]
[
0
0
0
]
T
Desired attitude angle[φdθdψd]
[
30
°
0
0
]
T
Initial angular velocityω0
[
1.1
1.2
1.6
]
T
°
/
s
Desired angular velocityωd
[
0
0
0
]
T
°
/
s
Initial gimbal angleδ0
[
0
0
0
-
90
°
]
T
Maximal angular velocityωm
8
°
/
s
Maximal control torqueum
1.5 N·m
Maximal gimbal angular velocityδ˙m
10
°
/
sThe simulation results are shown in Figures3, 4, 5, 6, 7, and 8. The solid line denotes the desired path planned by the improved GA. The dotted line expresses the results of the practical numerical simulation. The planning time and the whole maneuvering time are 7 s and 12 s, respectively.Euler angle of agile satellite.
(a)
(b)
(c)Angular velocity of agile satellite.
(a)
(b)
(c)Control torque.
(a)
(b)
(c)Gimbal angle of DGCMGs.
(a)
(b)
(c)
(d)Gimbal angle velocity of DGCMGs.
(a)
(b)
(c)
(d)Figure 8
Singularity measurement.The solid line in Figures3–8 indicates that the proposed GA is effective. The trajectory planned by the improved GA not only satisfies the desired maneuver ability, but also meets various physical constrains, such as input saturation, actuator saturation, angular velocity limit, and singularity measurement limit.Figures3–5 shows that the practical maneuvering curves of the state variables of the agile satellite are essentially consistent with the planning path. The results indicate good tracking performance of the derived controller.The conventional controller without path planning will cause the actuator saturation at the initial time because the initial error is large in the rapid maneuvering process. However, the results in Figure5 show that, the proposed control law not only keep the control torque away from the input restriction um, but also guarantee its smoothness.From the results in Figures6–8, there is few deviation of state variables curve of DGCMGs in the practical maneuvering process from the desired path. These errors are caused by using the attitude quaternion and velocity as reference trajectory. However, Figure 8 indicates that the singularity measurement curve is smooth and always away from the singularity in the practical maneuvering process. It is carried out to verify the ability to avoid DGCMG singularity of the proposed control law.


## 5. Conclusion

Agile satellite requires large-angle and rapid attitude maneuver ability. CMGs are expected to be applied to attitude control actuators of agile satellite. However, in the existing satellite attitude control systems using CMGs, the singularity problem cannot be considered from a global perspective. To solve the problem, an attitude tracking control algorithm with path planning is provided. Firstly, an improved GA is proposed to plan an energy optimal attitude path. Moreover, the desired path planned by GA satisfies certain maneuverability and various physical constrains, including input saturation, actuator saturation, angular velocity limit and singularity measurement limit. Then an attitude tracking controller is designed based on the sliding mode control method and the backstepping control method using hyperbolic tangent function as the nonlinear tracking function. And the adaptive control method is utilized to estimate uncertainties that results from the external disturbance and the inertial properties. From the results of the numerical simulations, the robustness and good tracking performance of the derived controller as well as its ability to avoid DGCMG singularity are verified.



---

*Source: 342043-2012-10-24.xml*