Fault detection and isolation of control moment gyros for satellite attitude control subsystem

Mechanical Systems and Signal Processing 135 (2020) 106419

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Fault detection and isolation of control moment gyros for satellite attitude control subsystem Afshin Rahimi a,⇑, Krishna Dev Kumar b, Hekmat Alighanbari b a b

Dept. of Mechanical, Automotive and Materials Engineering, Faculty of Engineering, University of Windsor, Windsor N9B 3P4, Canada Dept. of Aerospace Engineering, Faculty of Engineering and Architectural Science, Ryerson University, Toronto M5B 2K3, Canada

a r t i c l e

i n f o

Article history: Received 17 March 2019 Received in revised form 22 August 2019 Accepted 4 October 2019

Keywords: Fault detection Fault isolation Control moment gyro Satellite attitude control Unscented Kalman filter

a b s t r a c t Control moment gyros, as one of the most commonly used actuators onboard satellites, are prone to faults and failures. The ability to detect faults, isolate their location, and identify their severity can enhance mission success rate and reduce maintenance and damage costs extensively. Therefore, in this paper, a model-based fault detection, isolation, and identification scheme is proposed and evaluated. Firstly, an adaptive threshold fault detection algorithm is proposed, using Unscented Kalman filters in conjunction with residual and innovation sequences. Secondly, a fault isolation and identification algorithm is proposed using a binary grid search method to adapt joint estimation Kalman filter’s covariance matrix once a fault is detected. Extensive Monte Carlo simulations are conducted to evaluate the performance of the proposed schemes for a 3-axis stabilized satellite with a four-single-gimbal control moment gyro cluster. Results show superior performance of the proposed methods with faster tracking and more accuracy. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction The most common actuator technologies for small satellites are reaction wheels (RW), thrusters, and magnetorquers, while control moment gyros (CMG) are found in larger satellites. A CMG is a momentum exchange device consisting of a flywheel that is gimballed about one, two, or three of its axes. Gyroscopic torques are generated as the angular momentum vector is rotated about axes perpendicular to the flywheel’s spin axis. Several researchers have examined the problem of fault detection, isolation, and identification (FDII) [1–7]. Gertler [1] has surveyed model-based methods for fault detection and concludes that major quality issues for failure detection algorithms are isolability, sensitivity, and robustness. Marzat et al. [2] have reviewed model-based fault diagnosis approaches for aerospace systems. More recently, Gao et al. [3,4] have comprehensively reviewed fault diagnosis approaches and their applications from the model and signal-based perspectives. Tidriri et al. [8] have also investigated features of different model-based and data-driven fault diagnosis and health monitoring individually as well as hybrid approaches that incorporate the advantages of each. The reviewed literature suggests advantages of using Kalman filters as small false alarm rate, short detection delay, robustness to model uncertainty, and isolation of simultaneous faults with the only limitation being restrictive Gaussian noise assumption. In addition, parameter estimation using Kalman filters is considered as an effective approach for structural damage detection with the limitation of non-applicability for on-line identification due to time delays. An exten⇑ Corresponding author at: University of Windsor, Department of Mechanical, Automotive, and Materials Engineering, 401 Sunset Avenue, Windsor, ON N9B 3P4, Canada. E-mail address: [email protected] (A. Rahimi). https://doi.org/10.1016/j.ymssp.2019.106419 0888-3270/Ó 2019 Elsevier Ltd. All rights reserved.

2

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

sive review is published in a three-part series by Venkatasubramanian et al. [9–11]. Part I, II, and III of this series provide a review of quantitative model-based, qualitative model-based and history-based approaches, respectively. Parameter estimation methods [12–15], provide fault detection and isolation simultaneously since they provide the time, source, and size of the fault by estimating parameters of the system. To the authors’ knowledge, there are few studies explicitly published on the fault detection, isolation, and identification of the CMGs. Available studies on CMGs are mostly concerned with the fault-tolerant control (FTC) and the design of a controller rather than establishing an FDII framework. Zhang et al. [16] proposed a fault-tolerant attitude control system (ACS) for spacecraft with single gimbal CMG (SGCMG) under partial actuator failure and actuator saturation. Zhang et al. [16] mainly focus on the design of the controller and the FDII of the actuator (i.e. CMG) itself is not studied. Therefore, the main difference between the work presented in this paper is the approach in handling detection, isolation, and identification of the faults in the system. Mihai and Romulus Lungu [17] have proposed an adaptive neural network (NN)-based satellite ACS using dynamic inversion technique where the satellite is controlled through a cluster of variable speed CMGs (VSCMG). Their approach solves the FDII problem implicitly by designing an FTC where proportional-integral-derivative (PID) control is designed based on the linearized model of the system and NN updates the system model affected by unmodelled disturbances and system deteriorations over time. A neural network with at least one hidden layer is known to be a universal approximator, however, it comes at the price of extensive network training requirements. Therefore, one of the main disadvantages of using NN for system identification and FDII is the computational burden imposed on the system and the consequent delays that may arise from extensive training sequences. Zhang et al. [18] have proposed a faulttolerant decoupling control of spacecraft with SGCMGs based on active-disturbance-rejection control. The FTC employs two loops, one outer-loop to obtain the control torque where the actuator faults are not considered, and an inner-loop to calculate the virtual gimbal rate vector together with a singular direction-avoidance steering law. It is clear that the actuator faults are not estimated explicitly, and the controller is compensating for faults in the actuators if they arise. The issue with this approach from the FDII perspective is that it does not provide any information as to the location, size, and severity of the fault in the system. Therefore, to better capture the essence of the faults an explicit FDII is required. The few published literature where the FDII framework is explicitly proposed for CMGs includes Oh and Kim [19] where they have used a sliding mode observer (SMO) for detecting secular and oscillating disturbances. In their approach, a continuous second-order SMO is suggested to detect disturbances in the spacecraft ACS. The advantage of their approach is that both secular and oscillating external disturbances are aptly detected even in the existence of noise in gyro sensor data and the magnitude of disturbance is identified successively using frequency compensation of filtered gyro data. The internal wheel faults that arise from speed error are suggested to be obtained through their proposed method and applying a proper identification logic. The shortcomings with this work can be listed as twofold (1) A low-pass filter is used to filter out noise from gyro sensor data. A low-pass filter may not be sufficient to filter out noise, and other methods including Kalman filters can be employed for more robustness. (2) Fault isolation and identification cannot be explicitly achieved as it is merely an interpretation of reconstructed and filtered gyro sensor data where effect of multiple gyros on the satellite attitude can result in false alarms. Lee et al. [20] have extended the work in [19] by including CMG unit as hardware-in-the-loop in their system. Choi et al. [21] have used fuzzy logic and Q-learning for FDII of a small CMG-actuated satellite. The fuzzy logic is used to handle the nonlinear structures of the CMG while the Q-learning is employed for its online learning capabilities. One of the major drawbacks of the fuzzy logic is that the rules of combining membership functions (minmax rule) for conjunctive (AND) and disjunctive (OR) reasoning are not robust at all. Choi et al. have tried to compensate for this by employing Q-learning where Q-Learning is a fresher form of reinforcement learning algorithm that does not need a model of its environment and can be used online. However, the main problem of native Q-Learning is taking more time to reach the optimal value in normal circumstances with the obtained rewards [22]. Another drawback of their work is the adverse transient in switching the control strategy from normal to reconfigured after a fault is identified. In addition, developing the fuzzy logic for the system requires an extensive a priori knowledge to construct a full classification tree for motor faults as described in [21]. Shen et al. [23] have proposed an active FTC for spacecraft attitude maneuvers and actuator saturation and faults. For the fault estimation, they only estimate the total fault effects in place of individual faults. This means that the individual CMG faults are not captured and rather their collective effect on satellite attitude is estimated. The clear drawback of this approach from the FDII perspective is the inability to pinpoint explicit fault parameters in each actuator. In the published literature, there has not been much emphasis on the explicit FDII framework for CMGs as actuator systems onboard satellites. The main advantage of an explicit FDII is the ability to get more insight into the fault behavior. This information can be used to reconfigure system subcomponents (i.e., control logic, etc.) to adapt to the faulty condition. In addition, one can use this information to plan further developments on the faulty unit to avert such fault/failure in subsequent productions. Therefore, it is rational to develop an explicit FDII for CMGs that (1) Does not rely on the control system to compensate for the fault in actuator units [16–18]. (2) Employs more robust filtering techniques for reliable sensor readings, state/parameter estimation, and feature extraction [19,20]. (3) Does not require extensive training and re-training of the underlying models [21]. (4) Provides individual fault location, size, and severity for each actuator unit onboard the satellite rather than a lumped fault estimation [23]. Therefore, the main contributions of this work corresponding to the drawbacks reviewed above can be listed as (1) A new hierarchical model-based FDII approach is proposed for fault detection, isolation, and identification of the CMG units onboard satellite. (2) The proposed method does not require extensive a priori knowledge about the system for detection and isolation while (3) reduces the detection, isolation, and identification time delays and improves the detection and isolation accuracy compared to previously developed methods.

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

3

Therefore, the objective in this paper is to propose a new hierarchical FDII scheme based on adaptive thresholding detection and binary grid search covariance adaptation to remedy the shortcomings mentioned above. Furthermore, the performance of the proposed method needs to be evaluated on a 3-axis stabilized satellite in orbit using CMGs as its actuators. The contents of this paper are organized as follow. In Section 2, the problem definition is provided. In Section 3, the proposed methods are explained. In Section 4, the case study of 3-axis stabilized satellite with CMGs onboard as actuators is provided for performance evaluation. Section 5 lays out the simulation setup and provides detailed results and discussions on them. Section 6 provides concluding remarks, and finally, the Appendix includes the required background material for self-sufficiency of the paper. 2. Problem definition In this section, the problem of fault diagnosis in components of a general nonlinear system is formally stated. Consider the following nonlinear system described in discrete-time state-space representation:

8
 n > < nkþ1 ¼ f nk ; uk ; hk ; wk X : hkþ1 ¼ hk þ whk > : yk ¼ g ðnk ; hk Þ þ v k

ð1Þ

where nk 2 Rn is the state vector at time step k, uk 2 Rm is the control input vector, hk 2 Rl is the system parameter vector, yk 2 Rm is the measurement vector, wnk 2 Rn is the additive process noise for states, whk 2 Rl is the additive process noise for parameters, v k 2 Rm is the additive measurement noise. f ðÞ is a nonlinear process model, and g ð) is a nonlinear measurement model where in the case of full state measurement m ¼ n and yk ¼ nk þ v k is considered. The objective is to design and develop an FDII scheme that is capable of autonomously detecting the presence, isolating the location, and identifying the severity of faults in the system under the following assumptions: Assumption (i). The control signal and the state vector remain bounded prior to and after the occurrence of a fault. This is a reasonable assumption since the challenge in the FDII field is to detect low severity faults and anomalies since they do not drive the states or control signals to grow unbounded. Assumption (ii). There are no fault occurrences from the inception of a fault to its complete identification. This assumption ensures that during the FDII process, new changes in system performance do not shadow the state of the system under which the FDII module was activated and hence, result in misclassifications. This means that if a fault scenario is interrupted with another fault scenario (new faults occur before completion of the first FDII process) while the FDII is in the process of completing detection, isolation, and identification, the results from FDII may be corrupt, and misclassification can occur. Assumption (iii). Fault severity changes are ‘‘slow’’ compared to the system output dynamics. This assumption covers abrupt faults since once they occur, their severity does not change over time. In addition, incipient faults occur due to wear and tear, and their growth rates are often much slower than the system dynamics and hence, are covered under this assumption. Assumption (iv). All additive noises in the system are bounded and in the form of white Gaussian noise. The assumption on the white noise is made based on the classical definition of estimation with Kalman filter [24], however, in case of colored noise, augmented noise can be generated to fulfill the classical white noise requirement as discussed in [25]. 2.1. Formal representation of the system and fault models In this paper, it is assumed that the system component faults are reflected as changes in the physical system parameters [26]. Therefore, detecting changes in system parameters can be representative of faults in the system. The faulty system can be described the same as in Eq. (1) also known as multi-parameter fault model with:

hk ¼ h0 þ ak

ð2Þ

where h0 2 Rl is the nominal parameter values vector and ak 2 RL is the fault parameter vector containing L fault elements. The fault model given by Eq. (2) enables one to state the problem of nonlinear fault diagnosis in the form of an online nonlinear parameter estimation problem. 3. Methodology The proposed FDII in this study uses a hierarchical approach. The two levels of this hierarchy are (1) fault detection, (2) fault isolation and identification. Further detail about each level is provided in the following sections. The hierarchical approach ensures resources are allocated for computation only when required. Fig. 1 illustrates the overall hierarchy of the proposed methodology.where uk1 is the control input at the previous time step. yk is system output. nm k1 is the measured states vector in a single time step before fault detection and is used for fault reconstruction in the fault isolation and iden-

4

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

Fig. 1. Proposed hierarchical FDII framework.

tification module, Xs is the output of the FDII as the new system model, and z1 is the zero-order-old block that that at each time step provides the values from the previous time step. 3.1. Fault detection Fig. 2 shows the first level of the proposed methodology, which is fault detection. As can be seen, at this level, one model s ^m are the estimated states from the is running in parallel to the system where ysk is the output from the system. b n k and n k s unscented Kalman filter (UKF) based on Xs for system and model, respectively. K k and K m k are the calculated filter gains at time step k for the system and model, respectively. At the beginning of the simulations, it is assumed that Xs is the same as the system model and when faults occur in the system, the fault isolation, and identification module updates Xs to match ^m ^sk and y the current state of the system. y k are the estimated outputs for the system and model, respectively and r k is the generated residual at each time step. The function of the UKF at this level is to filter out noise and disturbances, eliminate unknown initial condition effect (for unmeasured states) and make the state estimation robust to some levels of modeling uncertainty for more reliability and accuracy in the estimated states. After r k is generated as

^k  y ^k rk ¼ y s

m

ð3Þ

scalar b is defined as [27]

b ¼ r Tk r k

ð4Þ

An adaptive threshold for anomaly detection (b0 ) [27] is used based on the available data and a fault is detected when b > b0 . After a fault is detected, the next step of the process, is to update Xs to match the current state of the system. This leads to the fault isolation and identification stage, which is detailed in the next section. 3.2. Filter adaptation algorithm The posterior estimates covariance matrix Pnn elements show available confidence in the estimated values for the system under study [28]. As the filter converges to the system parameters, the diagonal elements automatically decrease to show

Fig. 2. First level of the proposed hierarchical FDII framework.

5

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

more confidence in the estimation. However, when sudden changes in the system parameters occur, the filter cannot quickly track these changes due to increased confidence in its estimates. To remedy this problem, when a fault is detected, P nn elements need adjustment to compensate for the filter latency. In order to adjust Pnn elements, Rahimi et al. [27] proposed the algorithm in Table 1. where PP k is the temporary posterior estimate covariance matrix used for finding diagonal elements of the actual Pk matrix. The grid search limits for each PP k diagonal elements are based on numerical simulations where lower bound cannot be 0 due to numerical stability and upper bound needs to be large enough to speed up the tracking speed. The reason for choosing the maximum limit in the loop for 0.1 was through Monte Carlo analysis that found the higher values for the limit do not enhance the performance while extending the computation time for extra steps; therefore, 0.1 was chosen for the upper limit threshold. This method works well when the system is open-loop controlled and the system control input is pre-defined. However, when the system is closed-loop controlled, this approach does not yield satisfactory estimates. In addition, running an extensive grid search costs computation and the computation cost grows quadratically with the dimension of the system parameters that need to be estimated. This growth imposes numerical difficulties and delays in filter adaptation in case of faults. To remedy this problem, an enhanced version of this approach is proposed and tested in this paper to reduce the computation cost for the adaptation algorithm while maintaining its performance. The proposed method is referred to as binary grid covariance AUKF (GAUKF) shown in Table 2 where the idea is to only consider the binary choice of ‘‘increase” or ‘‘no change” for the diagonal elements of the posterior estimate covariance matrix. Where PPi is the binary combination to add to the previous step’s posterior covariance Pk1 , PP ki is the candid superimposed posteriori covariance, and rP i is the residual for estimates of the next time step using PPK i as posterior covariance. The main difference between CAUKF and GAUKF is the dimensionality reduction in GAUKF. The number of possibilities in CAUKF with the dimension of the parameters to be estimated as D and the number of grids choices for each element as n would be nD whereas for GAUKF this number would be equal to 2D . It is evident that there is a noticeable reduction in the curse of dimensionality in GAUKF. The rationale behind GAUKF still capturing the changes in system parameters and adapting filter to it is simple, whenever a change is detected in the fault detection step, it is in one of the parameters, by creating a grid of all possible options in step 2 of GAUKF we check increase in which parameters of the posterior covariance results in better estimates of the next time step or in other words, which elements have experienced change. Using this approach, we reduce the number of options needed to be examined in CAUKF by a notable amount while maintaining the performance of the estimator. And example of the process in GAUKF is shown in Eq. (5).

     1 0 1:01 0:003 0:01 0:003 þ ¼ 0 0 0:003 0:002 0:003 0:002 |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 

P k1

PPi

PPk

ð5Þ

i

In this example, which has the dimension of 2, for the combination added to the previous time step posterior covariance P k1 is assuming change only in the first estimated parameter and hence, only the element in row 1 and column 1 of PPi has the value of 1. In this case, because there are two parameters to be estimated, the total number of choices to loop through is 22 ¼ 4 cases only, while for the same scenario in CAUKF this number would equal to 102 ¼ 100. This further points out the computational advantage of using GAUKF over CAUKF. 3.3. Fault isolation and identification At this level, after a fault is detected, due to the nature of the proposed method, one can isolate and identify the fault in the system using Eq. (2). We can obtain the number of possibilities corresponding to the sum of all k-combinations of the faulty units using

L¼

n   X n p¼1

p

¼

n X p¼1

n! p!ðn  pÞ!

ð6Þ

where p is the number of faulty units at time step k and n is the number of available units that could become faulty. L determines the number of possibilities for the faulty system condition, each representing a combination of faulty units that could possibly be the current state of the system. If we run a dual sate/parameter estimation alongside the system after a fault is detected, we can estimate the fault parameters in the system and achieve the fault isolation and identification in one stage. To do this, first, we run a joint estimation to estimate system states and its parameters of interest. This helps to identify the severity of the fault from the available input/ outputs. Thus, the goal is to estimate parameters ak in Eq. (2) using a dual unscented Kalman filter estimation. Details of the dual estimation theory are explained in Appendix A. Fig. 3 illustrates the flow of this stage, where the main difference from the first stage in Fig. 2 is the addition of parameter estimates. Here, JUKF refers to the joint or dual UKF for joint estimation. s m As can be seen in Fig. 3, at this stage, system parameters are estimated along with its states. Therefore, b h and b h are k

k

introduced as parameter estimates for the system and model, respectively. Estimations are done for Xs . After parameters are estimated, a residual is calculated using normalized version of Eq. (4) as b ¼ r T Mr, where M is the weight matrix for nor-

6

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419 Table 1 Covariance adaptation algorithm (CAUKF) [27]. 1: 2: 3: 3a: 3b: 3c: 3d: 4: 5:

PP k ¼ P k1 and r old ¼ 1 For i = 1 to rows (P k1 ) For j = 0.1 to 1 step 0.1 SetPP kii ¼ j Calculate nk usingPP k Calculater new ¼ ðzk  nk ÞT ðzk  nk Þ If rnew < r old set P k ii ¼ j;r old ¼ r new Go to 3 Go to 2

Table 2 Binary grid covariance adaptation algorithm (GAUKF). 1: 2: 3: 3b: 3c: 4: 5:

rold ¼ 1 Create binary choices for PP k diagonal elements (0,1) For PPi ¼ 1tocombinations count (2D Þ PP ki ¼ P k1 þ PP i Calculate nk usingPP k Calculater P i ¼ ðzk  nk ÞT ðzk  nk Þ Go to 3 P k ¼ PP ki ðindexof ðminimumðrP i ÞÞ

Fig. 3. Third level of the proposed hierarchical FDII framework.

malization of different states in the system output yk and needs to be tuned for each system accordingly. This process is repeated until the residual is smaller than a set threshold of er ¼ 1  103 based on simulation results. When this criterion is met, the output of this stage would be the updated Xs model to replace the previous model in the first stage of the process and the monitoring continues. It should be noted that since the system parameters are directly estimated at this stage, one can inherently achieve fault isolation by evaluating the value for ak . In the nominal case, ak ¼ 0 and hk ¼ h0 . However, when a fault occurs, ak –0. Therefore, using a 1-to-k encoding, once can map binary faulty-healthy vectors and map them to the scenarios, as explained in Section 4.2.4. 4. Application case study: Fault diagnosis in CMG cluster of a satellite ACS To evaluate the performance of the proposed FDII, an ACS of a three-axis stabilized satellite using a four SGCMG (4SGCMG) cluster, is considered where actuator components experience faults. The simulation model consists of the nonlinear satellite attitude dynamics, a CMG cluster and a sliding mode controller (SMC) that is employed to stabilize the satellite. Fig. 4 depicts the block diagram of the closed-loop representation of the ACS subsystem that has been simulated in this paper. Following is the description of each module employed in this simulation: 4.1. Satellite dynamics The rotational equations of motion for a spacecraft actuation by momentum exchange devices such as CMGs can be described as [29] B H_ BI þ xBBI  HBBI ¼ se

ð7Þ

where xBBI is the angular velocity of the spacecraft relative to the inertial frame expressed in the body frame. se 2 R31 is the external torque. HBBI is the total angular momentum of the spacecraft relative to the inertial frame expressed in the body frame given by

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

7

Fig. 4. Proposed FDII simulation setup.

HBBI ¼ J xBBI þ h

ð8Þ T

where J is defined as J ¼ J s  AJw A and J s 2 R is the moment of inertia of the spacecraft including the actuators. J w 2 R44 ¼ diag ð½J w1 ; J w2 ; J w3 ; J w4 Þ denotes the axial moment of inertia of each momentum exchange device. h is the total momentum provided by the momentum exchange device expressed in the spacecraft body-fixed frame and A is the mapping matrix to map CMG’s actuator torques into spacecraft body axes. If CMGs are used as actuators, from Eqs. (43) and (44) we have [30] 33


 _ BBI ¼ xBBI  J s xBBI þ hCMG  h_ CMG þ se Jx

ð9Þ

where hCMG is defined in Eq. (61) and h_ CMG is defined in Eq. (62). 4.2. Attitude kinematics The kinematic equations for the spacecraft using quaternions can be formulated as



q_ v q_ 4



¼

  1 q4 I þ qv xBBL 2 qTv

ð10Þ

 

 T qv where q ¼ e sin U2 ; cos U2 ¼ is unit quaternion, U denotes the principal angle, e ¼ ½e1 ; e2 ; e3 T is the principal axis q4
 from Euler’s theorem e21 þ e22 þ e23 ¼ 1 . q4 2 R and qv 2 R31 ¼ ½q1 ; q2 ; q3 T denote the Euler parameters representing the spacecraft body frame orientation with respect to the orbital frame where qTv qv þ q4 ¼ 1. I 2 R33 is the identity matrix and q v is the skew-symmetric matrix of the quaternion vector (see Appendix B). 4.2.1. Actuators The selection of CMGs as actuators is well justified due to their popularity in the active satellite attitude control field. The CMG angular momentum is, in general, a function of CMG gimbal angles and flywheels angular speed where for RW only flywheel angular speeds are considered. The advantage of a gimbal angle addition is the added torque generation envelope for CMGs compared to RWs. The added gimbal requires an additional control component, which is referred to as the steering logic. The steering logic is responsible for converting required torque from the controller to the required gimbal angle rates that would generate that torque. In addition, the steering logic is responsible for avoiding singularities in the transformation matrices while requesting gimbal angle rates. The steering logic employed in this paper is adapted from Wie et al. [30]. For more detailed explanation on CMGs and steering logic see Appendix C. 4.2.2. Controller To obtain the desired attitude of qd 2 R41 and xd 2 R31 , a simplified version of a nonlinear sliding mode controller is adapted from [31] with the following quaternion tracking error

qe ¼ qd4 qv  q4 qdv þ qv qdv qe4 ¼ qd4 q4 þ qTdv qv where

qTe qe

Ce ¼

þ




q24e

q24e

ð11Þ

¼ 1. The corresponding rotation matrix C e ¼ C ðqe ; q4e Þ is given by

  qTe qe I þ 2qe qTe  2q4e qe

ð12Þ

8

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

31 where C Te C e ¼ 1, kC e k ¼ 1, det ðC e Þ ¼ 1, C_e ¼ x e C e , and I is the identity matrix. Next, the relative angular velocity xe 2 R is defined as follows

xe ¼ xBBL  C e xd

ð13Þ

Given these error terms, the sliding manifold is defined as

r ¼ xe þ kctr sgnðq4e Þqe

ð14Þ

where kctr > 0 is the sliding gain and sgnðq4e Þ is the sign function for q4e . The required control command is then obtained from

ur ¼ p0 r

ð15Þ

where p0 is a positive constant. With the required control command (ur ) and the steering logic, it is possible to calculate the gimbal angle rates as expressed in Eq. (65). In this work, all steering logic parameters are adapted from [30], kctr ¼ 1 based on values given in [31], p0 ¼ 0:1 based on the simulation results. 4.2.3. Fault formulation Four identical control moment gyros are used in a three-axis stabilized satellite in low earth orbit (LEO). The CMG cluster considered for this study is the pyramid configuration (see Appendix C, Fig. 10). Intermittent time-varying faults are injected into the mapping function of CMG (ACMG ) in Eq. (63) as [32]

ACMG ¼ ACMG F p

ð16Þ

with the CMG effectiveness matrix defined as

2

f p1 6 0 6 Fp ¼ 6 4 0 0

3

0

0

0

f p2

0

0 0

f p3 0

0 7 7 7 0 5

ð17Þ

f p4

where f pi corresponds to the effectiveness of the ith CMG and ranges from 0 to 1. Each fault parameter is an indicator of how well each CMG is performing with f p ¼ 1 for ideal performance and f p ¼ 0 for failure [33]. In other words, the multi-parametrized fault model is obtained by replacing f pi with f pi  ai , and ai are unknown fault parameters that indicate the possible presence of faults in each CMG unit. Due to the additive form of the fault parameters introduced above, the value for ai in the healthy condition would be zero and at any given time, the deviation from zero for any of these fault parameters could potentially be an indicator of the severity of the fault and its size. 4.2.4. Fault scenario To calculate the total number of combinations for a 4SGCMG assembly, we use Eq. (6) and it sums up to 16 unique combinations. To be able to refer to each combination (fault scenario) more conveniently, all possible cases are assigned with a number in Table 3 where the faulty unit number is defined in Fig. 10(a). While isolating a fault, the isolation and identification module would create a 1-to-k encoding of the faulty units and compare it with the Table 3 elements to output a faulty case. For example, if a1 –0, a3 –0 it means the corresponding CMGs are faulty and the resulting 1-to-k encoding vector would be [1 0 1 0]. This vector reflects scenario 6 in Table 3, hence, the isolation and identification module outputs scenario 6 as the current state of the system. The two stages of the proposed fault diagnosis in Fig. 1, are formulated for this case study as follows: 4.2.5. Fault detection In the first stage shown in Fig. 2, the satellite equations of motion are used as defined in Eqs. (9) and (10). The mapping matrix for the assembly configuration used in this study is defined in Eq. (63). The controller used to stabilize the attitude

Table 3 Fault scenarios for CMG cluster. Scenario no.

Faulty wheels

Scenario no.

Faulty wheels

0 1 2 3 4 5 6 7

None 1 2 3 4 1,2 1,3 1,4

8 9 10 11 12 13 14 15

2,3 2,4 3,4 1,2,3 1,2,4 1,3,4 2,3,4 1,2,3,4

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

9

uses Eq. (15) to calculate the required control input and translates it to the required gimbal angle rates using Eq. (65). At this stage, the state, output, and control vectors are defined as

n ¼ ½q1 ; q2 ; q3 ; x1 ; x2 ; x3 ; d1 ; d2 ; d3 ; d4  y ¼ ½q1 ; q2 ; q3 ; x1 ; x2 ; x3 ; d1 ; d2 ; d3 ; d4 u ¼



sx ; sy ; sz



ð18Þ

where qi and xi are the satellite’s quaternion and angular velocities, respectively. di are the gimbal angles for each CMG unit. Finally, ui are the commanded torque in the x; y; z directions and can be used to calculate h_ cmg in Eq. (64) and then plug it into Eq. (65) to obtain commanded gimbal angle rates d_ i as an input to each CMG unit. Together, these states form the reference

model of the system (i.e. Xs ). Using the reference model Xs , the current model of the system would initially be equal to the Xs where ai are equal to zero. The estimated output vectors b y s and b y m for the system and model, respectively use Xs for the UKF process in Eqs. (28)–(42) to obtain the states using the state estimation case (no parameter estimations) and are used to calculate the residual rk in Eq. (3). The residuals are used to form scalar b in Eq. (4) and are monitored with adaptive thresholding [27] for anomaly detection. Once an anomaly is detected, the next stage of the FDII scheme (fault isolation and m identification) is activated using the last healthy estimated states of the system (b n ) before the detection of the fault. k1

4.2.6. Fault isolation and identification In the second stage of the process shown in Fig. 2, the fault parameters (f pi ) are estimated to update Xs . At this stage, the system states, system parameters, control inputs, and measurements would be as follows

n ¼ ½q1 ; q2 ; q3 ; x1 ; x2 ; x3 ; d1 ; d2 ; d3 ; d4  h i h ¼ f p1 ; f p2 ; f p3 ; f p4

ð19Þ

y ¼ ½q1 ; q2 ; q3 ; x1 ; x2 ; x3 ; d1 ; d2 ; d3 ; d4 

u ¼ sx ; sy ; sz

The models used in the dual UKF process in Eqs. (28)–(42) to estimate the states n and parameters h are the same as Section 4.2.5. It is evident that when estimating h, a can inherently be calculated using

ab ¼ bh  h0

ð20Þ

b and b where a h are the estimates of the fault parameter and system parameter, respectively. When all ai are estimated, the final model Xs is output back to stage 1 of the process and replaces the model before fault detection. The monitoring continues, until another fault is detection where the proposed FDII process in Fig. 1 is followed. 5. Results and discussions To evaluate the performance of the proposed FDII scheme, numerical simulations were conducted. Numerical integration of the states is done in MATLAB using classical Runge-Kutta method (RK4) [34] with a step size of T s ¼ 0:01s and the total simulation time = 200 s. The machine used for the simulations is an IntelÒ CoreTM i7-4790 CPU with 3.60 GHz computing power and 16.00 GB RAM. Other simulation parameters are listed in Table 4. where the initial and desired conditions are calculated based on Euler angles to quaternion conversion via 3–21 (yaw





pitch-roll) sequence. The initial conditions for roll, pitch, and yaw are equal to 90 ; 10 ; 10 , respectively. The desired con





ditions for roll, pitch and yaw are equal to 0 ; 0 ; 0 , respectively. Process noise standard deviation (rQ ) and measurement Table 4 Simulation parameters [35]. Parameter

Value

UKF parameters ([a; j; k) Process noise covariance

[1,2,2] Q 0 ¼ r2Q Inn

Measurement noise covariance

R0 ¼ r2R Imm 2 0:015 0 4 0 0:017 0 0 2 1 0 60 1 1  105 6 40 0 0 0 [0.696, 0.707, [0,0,0] rad/s [0,0,0] [0,0,0] rad/s

Satellite moment of inertia (Js )

RWs moment of inertia (Jw )

Initial [q1 q2 q3  Initial [x1 x2 x3 ] Desired [q1 q2 q3  Desired [x1 x2 x3 ]

3 0 0 5[kg:m2  0:020 3 0 0 0 07 7[kg:m2  1 05 0 1 0.0]

10

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

noise standard deviation (rR ) are used in the process, and measurement noise matrices for the UKF and their values are set to

105 and 104 , respectively. Additive states noise with a standard deviation rn is added to the system outputs before measurement to produce noisy measurements and its value is defined in the following sections where the signal to noise ratio (SNR) is discussed. Given the measurement vectors in Eqs. (18) and (19), noise is added in the form of white Gaussian noise with zero mean and standard deviation rn in the form of

rn ¼ ½rni ; i ¼ 1;    ; M

ð21Þ

where M is the number of measurements in either Eqs. (18) or (19). The values of rni are determined in a way to set the same SNR for all measured vectors. Since each measured quantity is in a different order of magnitude, it is essential to set rni values in this way to ensure consistency for simulations and results. An out-of-phase abrupt fault scenario is described in Table 5, where each CMG’s effectiveness is reduced during an interval, and all CMGs’ are back to fully effective after the faulty period. Simulation results and detailed discussions on them are provided in the following section. 5.1. Results In this section, first, the system response results are demonstrated, followed by each stage of the FDII under different SNRs and finally a comprehensive Monte Carlo (MC) simulation. 5.1.1. System response Fig. 5 illustrates the response of the system under a faulty scenario described in Table 5 for 45 dB of SNR. Fig. 5(a) and (b) show that the system becomes stable before 20 sec into the simulation. After 20 sec mark, at each timestamp for the faults described in Table 5, there is a temporary disturbance in the system response until fault is detected, isolated, and identified and finally the controller can compensate for the loss of effectiveness with the updated system model. Fig. 5(c) illustrates the variations in the CMG units’ gimbal angels, which pertains to the requested torque from the controller. As can be observed, after each fault, there is a sudden change in the gimbals commanded angles to compensate for the loss of effectiveness. The large variance observed between 120 and 140 sec for gimbal angles, quaternions, and angular velocities is due to the error accumulation with the onset of 100 sec where f p4 value is changed. Although the parameter can be quickly identified, the error accumulation in the controller can still result in undesired deviations from the nominal reference points. System’s responses for the same scenario under 55 dB and 75 dB of SNR were also obtained; however, due to lack of space, only worst-case scenario (45 dB) is presented here. The rationale for choosing these three noise levels is explained in [36]. 5.1.2. Fault detection Fig. 6 illustrates an example of how the fault detection works based on adaptive thresholding for the scenario in Table 5. As can be seen, the adaptive threshold identifies when a fault is detected through an amplified residual metric b. It is important to note that due to the design of the proposed FDII procedure, faults need to be out-of-phase for the algorithm to properly capture and accommodate them by updating the reference model Xs . Otherwise, since the reference model parameters are not aligned with the current state of the system, fault detection module may properly work, but fault isolation and identification modules would be at fault. 5.1.3. Fault isolation and identification Table 6 shows the confusion matrix for Monte Carlo simulation under the 45 dB of SNR. Simulation outcomes for the same scenario described above under 55 dB and 75 dB of SNR were obtained; however, only the worst-case scenarios are included here. In these tables, symbol ‘‘A?” stands for actual model referring to Table 3 and symbol ‘‘I;” stands for the identified model from Table 3. A confusion matrix (as shown in Table 6) provides a visual-numerical representation for the performance of a classification method on a particular problem. This table shows how well the method classifies a category correctly. If the matrix is a diagonal matrix with off-diagonal elements all equal to zero, the performance is perfect, however, with the decrease in values of the diagonal elements and increase in the values of the off-diagonal elements, the performance decreases. It can be

Table 5 CMG fault parameters in-phase abrupt fault profile. Time (s)

f p1

f p2

f p3

f p4

t < 20 20  t < 40 40  t < 70 70  t < 100 100  t < 150 t  150

1 0.5 1 1 1 1

1 1 0.1 1 1 1

1 1 1 0.1 1 1

1 1 1 1 0.1 1

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

11

Fig. 5. Attitude response of the system under 45 dB noise: (a) satellite quaternions (b) satellite angular velocities (c) CMGs gimbal angles.

seen from Table 6 that the diagonal elements have values around ~90, which means about 90% of the time, the method correctly classified the category or identified the faulty units correctly. The average accuracy for the confusion matrix was calculated using the following formula:

accuracy% ¼

traceðCMÞ NCM

ð22Þ

where the traceðÞ function sums the diagonal elements of the matrix, subscript CM refers to confusion matrix, and nCM is the number of diagonal elements in the confusion matrix. Since there are 100 trials for each case, the accuracy is automatically

12

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

Fig. 6. Trend of b and b0 for the scenario described in Table 5.

Table 6 Confusion matrix for CMG using GAUKF – 45 dB Noise. A? I;

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

88 1 1 1 1 1 1 1 0 1 1 1 0 1 1

0 92 0 1 0 1 0 1 0 1 0 1 1 1 1

1 0 91 1 1 1 1 0 1 0 0 0 1 1 1

1 0 1 87 1 1 1 1 1 1 1 1 1 1 1

1 0 1 1 89 1 1 1 1 0 1 0 1 1 1

1 1 0 0 1 91 0 0 1 1 0 1 1 1 1

0 1 1 1 1 1 90 1 0 0 1 1 1 0 1

1 1 0 1 1 0 1 92 1 0 1 1 0 0 0

1 0 0 1 0 1 0 1 91 1 0 1 1 1 1

1 1 1 1 1 1 1 0 1 88 1 0 1 1 1

1 1 0 1 1 0 1 1 1 1 88 1 1 1 1

1 1 1 0 1 1 1 1 1 1 1 88 1 0 1

1 0 0 0 1 1 1 1 0 1 0 1 91 1 1

1 0 1 0 1 1 1 1 1 1 1 1 1 89 0

1 1 1 1 0 1 1 1 0 1 1 1 0 1 89

calculated as a percentage. With this definition for accuracy, the results shown in Table 7 can be obtained for all confusion matrices in this study. It can be seen from Table 7 that the accuracy for isolation ranges from 83.8% to 97.13%. It is important to note that the CAUKF has a meager accuracy percentage, and that is due to its lower performance while employed in a closed-loop controlled system. In the simulations it was observed that when CAUKF is put in closed-loop system, due to the contributions of the controller and stabilized system regardless of the fault, the performance of the CAUKF was deteriorated. Table 8 lists the average execution time for the Monte Carlo simulation for each method. It can be seen CAUKF takes the most time because it has to loop through extensive grids, GAUKF takes less time than CAUKF because it is only doing a binary grid search; however, it takes more time than the UKF. It should be noted that the 6.7% increase (ð8:74  8:19Þ=8:19  100) in execution time for GAUKF compared to UKF comes with the advantage of 8.2% increase (ð93:02  85:98Þ=85:98  100) in the

Table 7 Comparison of fault isolation accuracy under different noise levels. Noise level (SNR)

UKF%

CAUKF %

GAUKF %

75 dB 55 dB 45 dB Average

88.73 85.40 83.8 85.98

77.13 72.33 79.60 75.27

97.13 92.33 89.60 93.02

Table 8 Comparison of execution time. Method

UKF

CAUKF

GAUKF

Execution Time (sec)

8.19

9.69

8.74

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

Fig. 7. Fault identification performance under 45 dB noise.

13

14

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

average isolation accuracy as shown in Table 7 for 45 dB. Therefore, the benefit of incorporating the proposed GAUKF is the improved accuracy and reduced computation time compared to CAUKF. For fault identification, fault parameter changes for each CMG unit are examined. The parameter estimation used here is dual estimation described in Section 4.2.6 [27]. Fig. 7 shows the results of the fault identification for the scenario described in Table 5, under 45 dB of SNR. It can be seen that the GAUKF estimates are the closest to the true values. In addition, large spikes in the CAUKF as well as complete diversions in UKF (for f p4 after 150 sec) are the visual indicators of the superior performance of the proposed GAUKF. Fig. 7 demonstrates that the effectiveness parameters (f pi ) are properly estimated. Furthermore, it shows how one can calculate the respective fault parameters ai are as the difference between nominal (i.e., zero) and current values. Using fault parameters values, one can determine the severity of the fault and its size. Performance of the fault identification module for the same scenario in Table 5 under 55 dB and 75 dB of SNR were also obtained; however, only the worst-case scenario is presented here. In Fig. 7 for f p4 , it can be seen that the estimates for f p4 after 150 sec in the simulation are not as accurate as the previous time steps. It is important to note that UKF completely fails to readjust to the new value, CAUKF diverges over time while GAUKF stabilizes with a bias. The inferior estimation performance for the f p4 compared to other f pi s is due to the lower sensitivity of the system dynamics to the outputs of fourth CMG. The estimates work based on the residual of the system compared to model, and when the residual is very small, the effect results in inferior estimates. Hence, the trend observed in Fig. 7 for f p4 . Root mean squared error (RMSE) is calculated using

RMSE ¼

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u tP

2 u end b u f  f pi;t tt¼tstart pi;t t end  t start

ð23Þ

where b f pi is the estimated parameter value, f pi is the true parameter value, n is the number of parameters, and t start and tend are the start and end simulation time for the RMSE calculation, respectively. The RSMEs over the pre-fault, faulty, and post-fault periods (corresponding to tstart and t end ) for the fault scenario in Table 5 are summarized in Table 9 under 75 dB, 55 dB, and 45 dB SNR levels. The values in Table 9 are calculated for each parameter f pi and averaged over the range before fault occurrence (pre-fault), during faulty conditions (faulty), and after the fault occurrence is recovered (post-fault). It should be noted that the fault period for each f pi is different and non-overlapping. However, the values presented in Table 9 are the average of all for each noise level. The results in Table 9 show the superior performance of the GAUKF compared to CAUKF and UKF. In addition, it is evident that the increase in noise level reduces the performance of the filter and hence increases the RMSE values. 5.1.4. Monte Carlo analysis A confusion matrix is formed under different fault scenarios for different SNR levels. The tables are filled with the data of 100 simulations for each fault scenario in Table 5 under random initial conditions for satellite attitude (q; xÞ, fault start time (TF s ), fault duration (TF d ) that determines how long the fault would remain since TF s , and fault magnitude (ai ). The range for each of these parameters is given in Table 10. These ranges are chosen to cover a reasonable array of possible conditions given the satellite CMG effectiveness. It should be noted that the fault start and duration times are chosen within the given range in such a way that they remain out-of-phase. This check is done after a value is chosen for the MC simulation and if the condition is not satisfied, the start time and duration of the fault are re-initiated until the out-of-phase is satisfied for all effectiveness parameters. The fault detection rate (FDR) is defined as true positives (TP), the number of correctly identified faults, divided by all cases. The false alarm rate (FAR) is defined as false positives (FP), the number of faults detected while the system did not experience any changes in its fault parameters, divided by all cases. The missed fault rate (MFR) is defined as false negatives (FN), the number faults undetected while the system experienced changes in its fault parameters, divided by all cases. With these definitions, Table 11 provides a summary of fault detection performance under different noise levels. Furthermore, correlations between the size of the fault and its corresponding fault detection times are illustrated in Fig. 8 under 45 dB of SNR. Fault detection time is defined as the time between when the change in system parameters occurs to the time that the fault detection module activates the fault isolation module after b passes b0 threshold. Fault isolation time is defined as the time between when the fault isolation and identification is activated until it outputs the corrected model Xs . Different SNR levels for these correlations are compared in Table 12. It is important to note that the GAUKF takes the same amount of time to isolate and identify faults since it has a pre-set binary combination to check for each fault. The CAUKF, however, needs to check a larger grid and compare each one to find the best estimate for the next step and hence, the variable minimum and maximum isolation and identification delays. 5.2. Discussions Based on the results provided above, the performance of the proposed FDII method is discussed under different influencing factors as follows:

15

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419 Table 9 Comparison of the f pi RSME under different noise levels. Method

Noise level (SNR)

Pre-fault (102 )

Faulty (102 )

Post-fault (102 )

UKF

75 dB 55 dB 45 dB

9.23 10.61 11.41

12.30 13.87 15.34

9.45 10.32 11.31

CAUKF

75 dB 55 dB 45 dB

3.12 3.54 3.87

4.45 4.87 5.45

3.24 3.76 3.41

GAUKF

75 dB 55 dB 45 dB

1.74 1.98 2.10

2.15 2.44 3.14

1.83 2.06 2.16

Table 10 Parameters range for Monte Carlo simulations.

Min Max

qi

xi (rad/s)

TF s (s)

TF d (s)

ai

1 1

0.5 0.5

0 199

0 199-TF s

0 1

Table 11 Comparison of the fault detection performance under different noise levels for GAUKF. Noise level (SNR)

TP

FP

FN

FDR (%)

FAR (%)

MFR (%)

75 dB 55 dB 45 dB

1500 1495 1489

0 2 5

0 5 11

100 99.67 99.27

0 0.13 0.33

0 0.33 0.73

Fig. 8. Effect of fault size on FDII time delays under 45 dB noise for fault detection for GAUKF.

5.2.1. Fault severity The scale at which the change in CMG effectiveness parameters affects FDII performance was evaluated. Fig. 8 illustrates the correlation between the fault size and detection, isolation, and identification time delays. These plots are discontinuous after 95%. This means that the effectiveness parameter (f pi Þ values cannot be smaller than 5% of their nominal values (1); otherwise, the FDII module may provide wrong results. It can also be observed that the larger the fault parameters (ai ),

Table 12 Comparison of the time delays in FDII modules under different noise levels. Noise level (SNR)

Limit

Detection (s)

Isolation and identification (s) UKF

CAUKF

GAUKF

75 dB

Min Max

0.04 0.21

0.01 0.01

0.02 0.36

0.02 0.02

55 dB

Min Max

0.05 0.31

0.02 0.02

0.03 0.49

0.03 0.03

45 dB

Min Max

0.06 0.47

0.04 0.04

0.05 0.53

0.05 0.05

16

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

the more time it takes for the FDII to detect, isolate and identify them. It should be noted that small changes in fault parameters (< 5%) will not have drastic effects on satellite performance or cause a system failure. As such, it can be argued even in the presence of small changes in RW parameters, the system would still try to complete its mission objectives, and hence, this shortcoming for the proposed FDII module would not put the system at risk. 5.2.2. Noise Effect of noise on the performance of the proposed FDII is evaluated under three different noise levels (75 dB, 55 dB, 45 dB). Attitude response for the system under these three noise levels shows that although noise is adding jitter to the system response, it is not affecting FDII and controller effectiveness. Furthermore, the performance of the fault identification module in estimating fault parameters under different noise levels shows that although a subtle jitter is added in the estimation results, the overall performance of the fault identification is robust to measurement noise. Table 9 summarizes RMSE for fault identification estimations under various noise levels for the scenario described in Table 5. It can be seen that RSMEs are within acceptable boundaries and that during the fault there are larger RSMEs, which increase with the noise level. Finally, the confusion matrices for different noise levels are summarized in Table 11 into performance metrics FDR, FAR, and MFR. It can be seen that FDR decreases with the increase in noise level and both FAR and MFR increase with noise. This is expected since noise introduces uncertainties in the system; however, the changes in performance are within a negligible bound (<1%). 5.2.3. Time delays Table 12 summarizes the time delays under different modules of the proposed FDII for different noise levels. The majority of the time is spent in the isolation module due to its iterative structure. Furthermore, it can be observed that an increase in noise level, increases the time delays, which is expected since more noise introduces more uncertainty and for each module of the FDII this can mean more time to compensate for that. One can see the difference between the minimum and maximum CAUKF and GAUKF by about by ~0.5 sec. This is due to the iterative process for the CAUKF where it has to loop through all grid options for the posterior estimate covariance matrix, while in GAUKF, the binary grid approach eliminates the need to go through each option and reduces it to only two options for each element. The UKF has the least delay among the three approaches since it is not employing any iterative processes in its estimate covariance matrix update step. Although this can lead to less delay in estimations, divergence can happen for sudden changes in the estimated value once the estimation has settled as discussed in the methodology section. 6. Conclusions In this paper, an adaptive binary grid scheme was proposed to adapt states/parameters covariance matrix of the UKF in the case of sudden changes in non-measurable system parameters. The adaptive scheme was developed based on innovation and residual sequences of system states. The proposed scheme was tested on a 4SGCMG configuration onboard a closed-loop controlled satellite, to evaluate its performance. The proposed scheme provides a means to implement integrated FDII where fault detection, isolation, and identification module works collectively with the other subsystems (i.e., control logic, etc.) to not only help avert faulty units ending in failures, but also provide informative data on faults location, size, and severity. Extensive Monte Carlo simulations comprised out-of-phase abrupt faults in all four CMGs where study of concurrent faults is set as future work. Faults were investigated with different inception times and different severity for each CMG, and confusion matrices were generated. Results showed close to 93% success rate in isolation results compared to 75% for CAUKF and 85% for UKF. It was also observed that the computation time was decreased, and performance was increased in GAUKF compared to CAUKF. It can be concluded that the proposed GAUKF scheme performed superior compared to UKF and CAUKF in terms of speed and accuracy while employed in a closed-loop controlled system. GAUKF can be implemented where a priori knowledge about system performance is not available and computational resources can allow for the additional calculation requirements. Therefore, small satellites with limited space for hardware redundancy could be a reasonable target to implement such method and develop fail-safe satellites. Acknowledgment The authors acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) and Vanier scholarship. Appendix In this section, the additional preliminaries required for understanding the content of this paper are presented. A. State and parameter estimation with UKF

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

17

State and parameter estimation with UKF under the assumption of full state measurement is described in this section. Consider the nonlinear system in Eq. (1) where the process and measurement noise are assumed to be uncorrelated zeroh i h i mean Gaussian white noise with covariances E wk wTj ¼ Q k dkj , E v k v Tj ¼ Rk dkj , respectively. The process noise covariance Q k is non-negative definite, the measurements noise covariance Rk is positive definite, and dkj is the Kronecker-d function. When the system equation is in continuous form as in Eq. (1), it can be transformed to the discrete domain using

Z

ðkþ1ÞT s

nkþ1 ¼ nk þ

f ðnðt Þ; uðtÞ; hðt Þ; wðtÞÞdt

ð24Þ

kT s

where T s is the discretization step size. The integration can be achieved using the fourth-order Runge-Kutta (RK4) method. To jointly estimate the state and parameters of the system the following augmented vector is defined:

h a T b xk ¼ b xk

iT

ð25Þ

wk T

h iT T T b where superscript T is used for vector transpose operation, b x k 2 Rnþl ¼ b has the state and parameter estimates, nk hk h iT respectively. The vector wk 2 Rnþl ¼ wnk T whk T contains the process noises in the evolution of n and h. Analogously, the augmented covariance matrix is defined as

"

Pak ¼

Pxk

Px;w k

Pkw;x

Pw k

#

ð26Þ

where P xk consists of the state and parameter error covariances, while Pw k includes the process noise covariance associated with the state and parameters. In addition, the off-diagonal entries of the cross-covariance terms denoted by P};} k . It is evident a that all elements of b x k and Pak are of appropriate dimensions. Given the initial conditions,

"

a b x0 ¼

# " b Px0 x0 a ; P0 ¼ 0 0w

0

#

ð27Þ

Pw 0

with the dimension of the augmented system state b x a equal to L ¼ n þ l, the following steps are taken to calculate the estimates and covariances in each iteration of the UKF process: 1. Compute weights ðmÞ

W0

¼

j

ðc Þ

Lþj

; W0 ¼

j

Lþj


 ðm Þ ðc Þ þ 1  a2 þ k ; W i ¼ W i ¼

1 2ðL þ jÞ

where W ðmÞ is a component weight for mean calculations, W ðcÞ is a component weight for covariance calculation, 0  a < 1, and k  0 are control factors for the spread of sigma points.

ð28Þ

j  0,

2. Establish symmetric sigma points about the state estimate




X ak1

 0

ba ¼X k1

a bj ¼ b X x k1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ffi ðL þ jÞ Pak1 j ; 8j ¼ 1; 2;    ; n

ð29Þ ð30Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
 where the ðL þ jÞ P ak1 j terms denote the scaled jth rows or columns of the square root of the augmented covariance matrix at the previous time step. after calculating sigma points, regroup them in the following matrix of L rows and 2L þ 1 columns



X ak1 ¼

X xk1



Xw k1

ð31Þ

were vxk1 contains the sigma point rows associated with states and parameters, and X w k1 contains the sigma point rows associated with the state and parameters process noises. 3. Propagate the sigma points through the nonlinear dynamic



X kjk1 ¼ F X ak1 ; uk1

ð32Þ

18

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

xk ¼

2L X

ðmÞ

Wi

ð33Þ

vi;kjk1

i¼0

where F ½ is the nonlinear dynamic model, iis the appropriate sigma point column index, and superscript ðÞ indicates the pre-process variable. 4. Compute the predicted covariance

Pk ¼

2L X

ðc Þ

Wi



X i;kjk1  b xk

 T

xk X i;kjk1  b

ð34Þ

i¼0

where P k is the augmented states and parameters covariance matrix. 5. Compute new sigma points Given the predicted mean in Eq. (33) and covariance in Eq. (34) recompute a new set of sigma points as defined in Eqs. (29)–(31)

h X k ¼ X xk T

Xw k

T

iT

ð35Þ

6. Instantiate the new sigma points through observation model



Y kjk1 ¼ G X k  b yk ¼

2L X

ð36Þ

ðm Þ

W i Y i;kjk1

ð37Þ

i¼0

where Y kjk1 is the result of instantiated sigma points through the observation model, G½ is the nonlinear observation model,  and b y k is the weighted mean of variable y. 7. Obtain the innovation and cross-covariance matrices

Pyk yk ¼

2L X

ðc Þ

Wi

  T Y i;kjk1  b y k Y i;kjk1  b y k þ Pv

ð38Þ

i¼0

Pyk xk ¼

2L X

ðc Þ

Wi



X i;kjk1  b xk

 T

yk Y i;kjk1  b

ð39Þ

i¼0

where Pyk yk is the innovation covariance matrix, Pyk xk is cross-covariance matrix, and P v is the measurement noise covariance. 8. Calculate gain and update

Kk ¼ Pyk xk P1 yk yk

ð40Þ


  b x k þ Kk y k  b yk xk ¼ b

ð41Þ

Pxk ¼ Pk  Kk Pyk yk Kk T

ð42Þ

where Kk is the Kalman filter gain, b x k is the new augmented states and parameters estimate, and Pxk is the updated states and parameters covariance matrix. B. Satellite attitude dynamics and kinematics The satellite is modeled as a rigid body. The coordinate frames used for satellite modeling are shown in Fig. 9. An Earthcentered inertial (ECI) reference frame denoted by I  X I Y I Z I has its origin at the center of the Earth, where Z I axis passes through the celestial North Pole, X I axis points towards the vernal equinox, and Y I axis completes the triad with the right-handed rule. The orbital frame also known as local vertical local horizontal (LVLH) is denoted by L  X O Y O Z O and is fixed at the center of the spacecraft. In LVLH, X O axis is in the direction of motion, Y O axis in the opposite direction of the

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

19

Fig. 9. Geometry of orbit motion for rigid body spacecraft.

angular velocity, and Z O axis is pointing at the Earth. True anomaly (for eccentric orbits) or reference angle h (for circular orbits) is measured from the line passing through the center of the Earth to the spacecraft. The body-fixed coordinate frame denote by B  X B Y B Z B has its origin at the spacecraft center of mass [37]. A vector of frame L relative to frame I, expressed in frame B is denoted by v BLI [38]. 9. Satellite dynamics The rotational equations of motion for a spacecraft actuation by momentum exchange devices such as CMGs can be described as [29] B H_ BI þ xBBI  HBBI ¼ se

ð43Þ

where xBBI is the angular velocity of the spacecraft relative to the inertial frame expressed in the body frame. se 2 R31 is the external torque. HBBI is the total angular momentum of the spacecraft relative to the inertial frame expressed in the body frame given by

HBBI ¼ JxBBI þ h

ð44Þ T

where J is defined as J ¼ J s  AJw A where J s 2 R is the moment of inertia of the spacecraft including the actuators. J w 2 R44 ¼ diag ð½J w1 ; J w2 ; J w3 ; J w4 Þ denotes the axial moment of inertia of each momentum exchange device. h is the total momentum provided by the momentum exchange device expressed in the spacecraft body-fixed frame. If CMGs are used as actuators, from Eqs. (43) and (44) we have [30] 33


 _ BBI ¼ xBBI  Js xBBI þ hCMG  h_ CMG þ se Jx

ð45Þ

where hCMG is defined in Eq. (61) h_ cmg is defined in Eq. (62). 10. Attitude kinematics The kinematic equations for the spacecraft using quaternions can be formulated as



q_ v q_ 4



¼

  1 q4 I þ qv xBBL 2 qTv

ð46Þ

 

 T qv is unit quaternion, U denotes the principal angle, e ¼ ½e1 ; e2 ; e3 T is the principal axis where q ¼ e sin U2 ; cos U2 ¼ q4
 from Euler’s theorem e21 þ e22 þ e23 ¼ 1 . q4 2 R and qv 2 R31 ¼ ½q1 ; q2 ; q3 T denote the Euler parameters representing the spacecraft body frame orientation with respect to the orbital frame where qTv qv þ q4 ¼ 1. I 2 R33 is the identity matrix and q v is the skew-symmetric matrix of the quaternion vector given by

20

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

2

q3

0

6 qv ¼ 4 q3

q2

7 q1 5

0

q2

3

q1

ð47Þ

0

The angular velocity of the body frame with respect to the inertial frame can be written as

xBBI ¼ xBBL þ xBLI

ð48Þ

The transformation matrix C BL is used to describe the orientation of the spacecraft body frame B with respect to the orbital frame L. It can be formulated in terms of Euler parameters as follows


 CBL ¼ q24  qTv qv I þ 2qv qTv  2q4 qv The orbital angular velocity expressed in the body framex defined as

ð49Þ B OI

(i.e., the speed the orbit is rotating with the Earth) can be

xBLI ¼ CBL ½ 0 x0 0 T

ð50Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffi where x0 ¼ h_ is the orbital frame angular velocity. This value is equal to h_ ¼ le =R3c for circular orbits were le is the gravitational constant parameter of the Earth equal to 398600km3 =s2 and Rc is the spacecraft distance from the center of the Earth. This will be equal to Rc ¼ RE þ hs where RE is the average Earth radius equal to 6378km and h is the altitude of the spacecraft from the surface of the Earth. In this study hs ¼ 500km. In orbital satellite motion, it is desired to work with the spacecraft body angular velocity relative to the orbital frame L _ BBL . From theory we have [39] expressed in the body frame x

h i B  C_ L ¼  xBBL CBL

ð51Þ

Taking a time derivative of Eq. (48) and applying it to Eq. (51) we have

h



i

x_ BBI ¼ x_ BBL  xBBL xBLI

ð52Þ

_ BBL we have Since we are interested in x

h



i

x_ BBL ¼ x_ BBI þ xBBL xBLI

ð53Þ

The external torque is assumed to include gravity gradient torque (sgg ), solar radiation pressure (ssrp ), aerodynamic torque (saero ), and other disturbances (sd ). se ¼ sgg þ ssrp þ saero þ sd . These disturbances can be defined as [35,40,41]

sgg ¼ 3x20 c3 Js c3 ; c3 ¼ CBO ½ 0 0 1 T

3 2 2  105 ½1  2 sin ðx0 tÞ Tsx 7 6

7 6 6 7 3 7 1  10 ¼ 4 Tsy 5 ¼ Sf 6 ð x t Þ cos 0 7 6 5 4

Tsz 5  5  10 cosðx0 tÞ 2

ssrp

saero



sd ¼

3

3 1 þ sin ðx0 tÞ þ 12 sinð2x0 tÞ

6

7 6 7 2 1 7 ¼ 4 Tay 5 ¼ 1:36Af 6 4  5  10 4 þ 2sinðx0 tÞ þ 2 sinð2x0 tÞ 5

Taz  1 þ sinðx0 tÞ þ 12 sinð2x0 tÞ 2

ð54Þ

Tax

3

ð55Þ

2

2 3  sinðx0 tÞ 1 2 6 7 þ kxBBO k 4 cosðx0 tÞ 5 2 cosðx0 tÞ

ð56Þ

ð57Þ

where positive scaling factors Sf and Af are calculated based on the worst expected disturbance torques in low earth orbit (LEO) and are equal to Sf ¼ 1:70  106 , Af ¼ 1:0  1011 for RyeSat (pico-satellite) [35]. C. Control moment gyro In general, the categorization of CMGs is based on the gimballing arrangements and the configuration used for redundancy management and failure accommodation. The two basic types of CMGs are (1) single-gimbal control moment gyros

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419

21

(SGCMGs) and (2) double-gimbal control moment gyros (DGCMGs). The CMG angular momentum is, in general, a function of CMG gimbal angles d ¼ ðd1 ;    ; dn Þ and flywheels angular speed X ¼ ðX1 ;    ; Xn Þ

HCMG ¼ Hðd; XÞ

ð58Þ

One of the approaches to the CMG steering logic is to use the differential relationship between gimbal angles and the CMG momentum vector. For such a method, the derivation of h is obtained as

h_ CMG ¼ ACMG d_

ð59Þ

where ACMG ¼ ACMG ðdÞ 2 R

3n



ACMG ¼

@h @hi ¼ @d @di



is the Jacobin matrix is

ð60Þ

The CMG steering logic then becomes the inverse of ACMG d_ ¼ h_ CMG . If we consider a typical pyramid mounting configuration for four SGCMGs as shown in Fig. 10 where each face is inclined with an out-of-plane angle of b ¼ 54:73deg, which makes the momentum envelope nearly spherical [42]. For the conventional pyramid mount (Fig. 10), the total CMG angular momentum can be expressed in the spacecraft reference frame as

2

hCMG ¼

4 X

6 hi ðdi ; Xi Þ ¼ 4

3

cbsd1

cd2

cbsd3

cd1

cbsd2

cd3

T 7 cbsd4 5  h01 ðX1 Þh02 ðX2 Þh03 ðX3 Þh04 ðX4 Þ

sbsd1

sbsd2

sbsd3

sbsd4

i¼1

cd4

ð61Þ

where hi is the angular momentum of each CMG expressed in the spacecraft reference frame. di are the gimbal angles, Xi are flywheel angular speed, and h0i is the momentum magnitude for the ith CMG. The time derivative of the CMG angular momentum can be obtained as

h_ CMG ¼

4 X



h_ i ðdi ; Xi Þ ¼ h01 ðX1 Þh02 ðX2 Þh03 ðX3 Þh04 ðX4 Þ ACMG d_

ð62Þ

i¼1

where d is the gimbal angle vector and

2

cbcd1 6 ¼ 4 sd1 sbcd1

ACMG

cbcd3

cbcd2

sd3

7 cbcd4 5

sbcd2

sbcd3

sbcd4

For a known control torque

sd4

3

sd2

ð63Þ

sc , the CMG torque command h_ is chosen as

h_ CMG ¼ u ¼ sc  xBBI  hCMG

ð64Þ

_ given h0 ¼ h0 ¼ h0 ¼ h0 ¼ h0 is obtained as [30] And the gimbal rate command d, 1 2 3 4

d_ ¼



 þ 1 h_ CMG A h0 CMG

Fig. 10. Pyramid mounting configuration of four SGCMGs (a) schematics (b) free body diagram.

ð65Þ

22

A. Rahimi et al. / Mechanical Systems and Signal Processing 135 (2020) 106419



1 T T where Aþ , often referred to as the pseudoinverse steering logic and most CMG steering laws deterCMG ¼ ACMG ACMG ACMG mine the gimbal rate commands with its variations. References [1] J.J. Gertler, Survey of model-based failure detection and isolation in complex plants, IEEE Control Syst. Mag. 8 (6) (1988) 3–11. [2] J. Marzat, H. Piet-Lahanier, F. Damongeot, E. Walter, Model-based fault diagnosis for aerospace systems: a survey, Proc. Inst. Mech. Eng. Part G J. Aerosp. Eng. 226 (10) (2012) 1329–1360. [3] Z. Gao, C. Cecati, S.X. Ding, A survey of fault diagnosis and fault-tolerant techniques Part I: fault diagnosis with model-based and signal-based approaches, IEEE Trans. Ind. Electron. 62 (6) (2015) 3757–3767. [4] Z. Gao, C. Cecati, S. Ding, A survey of fault diagnosis and fault-tolerant techniques part II: fault diagnosis with knowledge-based and hybrid/active approaches, IEEE Trans. Ind. Electron. 62 (6) (2015) 3768–3774. [5] R. Isermannt, Process fault detection based on modeling and estimation methods – A survey, Automatica 20 (4) (1984) 387–404. [6] D. Kim, D. Choi, H.S. Oh, Failure diagnosis of actuator using the EKF and NP theorem, Int. J. Control. Autom. Syst. 11 (3) (2013) 450–459. [7] P.M. Frank, Fault diagnosis in dynamic systems using analytical and knowledge-based redundancy a survey and some new results*, Automatica 26 (3) (1990) 138–143. [8] K. Tidriri, N. Chatti, S. Verron, T. Tiplica, Bridging data-driven and model-based approaches for process fault diagnosis and health monitoring: a review of researches and future challenges, Annu. Rev. Control 42 (2016) 63–81. [9] V. Venkatasubramanian, R. Rengaswamy, K. Yin, S.N. Kavuri, A review of process fault detection and diagnosis: part I: quantitative model-based methods, Comput. Chem. Eng. 27 (3) (2003) 293–311. [10] V. Venkatasubramanian, R. Rengaswamy, S.N. Kavuri, A review of process fault detection and diagnosis: part II: qualitative models and search strategies, Comput. Chem. Eng. 27 (3) (2003) 313–326. [11] V. Venkatasubramanian, R. Rengaswamy, S.N. Kavuri, K. Yin, A review of process fault detection and diagnosis: part III: process history based methods, Comput. Chem. Eng. 27 (3) (2003) 327–346. [12] T. Jiang, K. Khorasani, S. Tafazoli, Parameter estimation-based fault detection, isolation and recovery for nonlinear satellite models, IEEE Trans. Control Syst. Technol. 16 (4) (2008) 799–808. [13] X. Li, Q. Zhang, H. Su, An adaptive observer for joint estimation of states and parameters in both state and output equations, Int. J. Adapt. Control Signal Process. 25 (9) (2011) 831–842. [14] A. Rahimi, K.D. Kumar, H. Alighanbari, Enhanced adaptive unscented kalman filter for reaction wheels, IEEE Trans. Aerosp. Electron. Syst. 51 (2) (2015) 1568–1575. [15] H. Ye, W. Wang, S. Zhai, Fault diagnosis based on parameter estimation in closed-loop systems, IET Control Theory Appl. 9 (7) (2015) 1146–1153. [16] F. Zhang, L. Jin, S. Xu, Fault tolerant attitude control for spacecraft with SGCMGs under actuator partial failure and actuator saturation, Acta Astronaut. 132 (2017) 303–311. [17] M. Lungu, R. Lungu, Adaptive neural network-based satellite attitude control by using the dynamic inversion technique and a VSCMG pyramidal cluster, Complexity 2019 (2019) 1–16. [18] F. Zhang, L. Jin, S. Xu, Y. Zhao, Fault-tolerant decoupling control for spacecraft with SGCMGs based on an active-disturbance rejection-control technique, J. Aerosp. Eng. 31 (2) (2018) 04018001. [19] H. Oh, J. Kim, ‘‘The detection and identification of control moment gyro faults,” in Second IAA Conference on Dynamics and Control of Space Systems 2014, DyCoSS’2014, 2015, pp. 1–10. [20] Jung-Hyung Lee, Hun-Jo Lee, Joon-Yong Lee, Oh. Hwa-Suk, External disturbance detection and its application on the identification of fault in CMG system, J. Energy Power Eng. 12 (6) (2018) 281–288. [21] Y.-C. Choi, J.-H. Son, H.-S. Ahn, Fault detection and isolation for a small CMG-based satellite: a fuzzy Q-learning approach, Aerosp. Sci. Technol. 47 (2015) 340–355. [22] M.S. Manju, An analysis of Q-learning algorithms with strategies of reward function, Int. J. Comput. Sci. Eng. 3 (2) (2011) 814–820. [23] Q. Shen, C. Yue, C.H. Goh, D. Wang, Active fault-tolerant control system design for spacecraft attitude maneuvers with actuator saturation and faults, IEEE Trans. Ind. Electron. 66 (5) (2019) 3763–3772. [24] R.E. Kalman, A new approach to linear filtering and prediction problems, J. Basic Eng. 82 (1) (1960) 35. http://FluidsEngineering.asmedigitalcollection. asme.org/article.aspx?articleid=1430402, https://doi.org/10.1115/1.3662552. [25] Y. Bulut, Applied Kalman Filter Theory, Northeastern University, 2011. [26] E. Sobhani-Tehrani, H.A. Talebi, K. Khorasani, Hybrid fault diagnosis of nonlinear systems using neural parameter estimators, Neural Netw. 50 (2014) 12–32. [27] A. Rahimi, K.D. Kumar, H. Alighanbari, Fault estimation of satellite reaction wheels using covariance based adaptive unscented Kalman filter, Acta Astronaut. 134 (2017) 159–169. [28] S.S. Bisht, M.P. Singh, An adaptive unscented kalman filter for tracking sudden stiffness changes, Mech. Syst. Signal Process. 49 (1–2) (2014) 181–195. [29] P. Tsiotras, H. Shen, C. Hall, Satellite attitude control and power tracking with energy/momentum wheels, J. Guid. Control. Dyn. 24 (1) (2001) 23–34. [30] B. Wie, D. Bailey, C. Heiberg, Singularity robust steering logic for redundant single-gimbal control moment gyros, J. Guid. Control. Dyn. 24 (5) (2001) 865–872. [31] K.D. Kumar, N. Abreu Godard, M. Sinha, Fault-tolerant attitude control of miniature satellites using reaction wheels, Acta Astronaut. 151 (May) (2018) 206–216. [32] N. Abreu, Fault Diagnosis with Adaptive Kalman Filters and CMG Design for Picosatellite ACS, Ryerson University, 2010. [33] D. Kim, J.D. Turner, J.L. Junkins, Optimal actuator failure control using a homotopy method, J. Guid. Control. Dyn. 38 (4) (2015) 623–630. [34] K.E. Arkinson, An Introduction to Numerical Analysis, Second ed., John Wiley & Sons, New York, NY, USA, 1989. [35] Godard, Krishna Dev Kumar, An-Min Zou, A novel single thruster control strategy for spacecraft attitude stabilization, Acta Astronaut. 86 (2013) 55– 67, https://doi.org/10.1016/j.actaastro.2012.12.018. [36] E. Sobhani-Tehrani, K. Khorasani, ‘‘Identification For Nonlinear Systems Using Hybrid Approach,” Master’s Thesis. pp. 12–13,15,18,37,92, 2008. [37] J. Li, K. D. Kumar, ‘‘Decentralized fuzzy fault tolerant control for multiple satellites attitude synchronization,” in 2011 IEEE International Conference on Fuzzy Systems (FUZZ-IEEE 2011), 2011, pp. 1836–1843. [38] D. Fragopoulos, M. Innocenti, Stability considerations in quaternion attitude control using discontinuous Lyapunov functions, IEE Proc. – Control Theory Appl. 151 (3) (2004) 253–258. [39] A.H.J. De Ruiter, C.J. Damaren, J.R. Forbes, Spacecraft Dynamics and Control An Introduction, first ed., John Wiley & Sons Ltd., 2013. [40] H. Schaub, J. L. Junkins, ‘‘Analytical mechanics of space systems,” AIAA Educ. Ser., p. 11,34,37,112,113, 2003. [41] W. Cai, X. Liao, D.Y. Song, Indirect robust adaptive fault -tolerant control for attitude tracking of spacecraft, J. Guid. Control. Dyn. 31 (5) (2008) 1456–1463. [42] C.H. Bong Wie David Bailey, Singularity robust steering logic for redundant single-gimbal control moment gyros, J. Guid Control 24 (5) (2001) 865–872.