{'fraction_non_alphanumeric': 0.024417507947043, 'fraction_numerical': 0.010130296573165195, 'mean_word_length': 4.582683307332293, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real-Time Tuning Unscented Kalman Filter for a Redundant Attitude Estimator in Microsatellites

This paper deals with microsatellite attitude determination systems which are a combination of a main estimator with high-power, high-precision sensors for higher accuracy estimation and a redundant estimator with low-power, lower-precision sensors for backup. Measurement data from all sensors in the redundant estimator are fused by the unscented Kalman ﬁlter to provide estimated attitude and the gyro bias values. Besides the accuracy of attitude sensors, the accuracy of this estimator depends largely on the selection of the process and measurement noise covariance matrices. In this paper, a novel real-time tuning unscented Kalman ﬁlter for redundant attitude estimator is introduced to tune these matrices eﬃciently in each ﬁlter step. The tuning process uses the estimated attitude of the main estimator as an independent truth reference data to calculate the cost function which is minimized by a downhill simplex algorithm. In the scheme developed in this paper a ﬁne-tuning process is used, which results in faster convergence speed and higher estimated accuracy of the redundant estimator. Another important feature of the developed ﬁlter is that a ﬂexibly estimated accuracy and system power consumption can be archived by choosing the duration and repeat frequency of turn-on time of the main estimator.


Introduction
Recently advanced microsatellites may use various types of attitude sensors in order to flexibly detect their attitude with the minimum power consumption. In the attitude determination system (ADS) of those satellites, a redundant estimator with low-power and lower-precision sensors is used as a backup estimator for a main estimator that uses high power, high-precision sensors for higher accuracy estimation. Although the main estimator with the highest estimated accuracy is always preferred. However, when the satellites work under power budget constraints, the redundant estimator should be used instead. Depending on the type of sensors which the main and redundant estimators use, the working area of these estimators may separate in terms of the balance of power and precision as in Fig. 1. Therefore, if the mission requirements for power consumption and estimated precision fall in the middle area, the ADS cannot archive these requirements. Under this point of view, the proposed filter is used to expand the working area of the redundant estimator and to provide a reasonably good estimator using additional data from the main estimator. The additional data can be discontinued in a sufficiently long time to reduce the system power consumption. The estimated attitude and gyro bias values are obtained by the filter algorithms. This estimated accuracy is improved by the use of a set of optimal parameters inside the filters which is updated in real-time throughout the tuning process. Before the detailed description of proposed method, it is necessary to give a brief review of attitude filters and parameter tuning problems.
The Kalman filter (KF) is a sequential state estimator that is optimal in the sense of minimizing the estimated error covariance. The selection of process noise and measurement noise covariance matrices Q and R is the primary factor that determines the steady-state value of Kalman gain K, 1) which decides the stability and performance of the filter. In the case of small gain K, measurements tend to be ignored and the system relies more on the model; otherwise, measurements have more effect on the performance than the model does. Therefore, Kalman gain K controls how the new measurement information is combined with the internal dynamic model and so determines the performance of the filter. 2) The original KF was designed for linear dynamical systems. In the case of nonlinear systems, such as the satellite attitude estimation problem, the KF was expanded to include nonlinear dynamics models, which led to the devel-opment of the extended Kalman filter (EKF) or unscented Kalman filter (UKF). 3) The EKF and UKF are both wellknown and flight-confirmed algorithms for satellite attitude estimation. [4][5][6] UKF has more advantages than EKF such as lower expected error and validity to higher-order expansions. Due to the linearization effect and the approximation of the state distribution effect in EKF and UKF, respectively, the determination of the Q and R matrices now becomes more difficult in practical implementation.
KF tuning refers to the process of selecting the elements of the Q and R matrices to improve the estimates with respect to a performance measure such as a cost function. It is difficult to determine the optimal values of all Q and R elements, especially in a real-time system, in which the noises of the model and measurement have a strong effect on the estimation. The conventional approaches to find the best parameter values of this tuning problem are trial and error or grid search methods. However, they require a lot of execution time and calculation efforts to work with multiple tuning parameters.
An automated method for KF tuning has been introduced by Powell. 2) The method formulates a stochastic cost function called performance index based on the true estimation error, typically obtained from extensive Monte Carlo simulations. The proposed method solves the associated minimization problem using the downhill simplex algorithm. The method mainly works for offline analysis due to the use of batch data and a large number of Monte Carlo runs.
Basil et al. introduced another tuning method working as an adaptive sequential estimation with unknown noise statistics. 7) In this method, the Q and R matrices are recalculated using N measurement points based on the minimization of the performance index that depends on the innovation vector which is the difference between the model-predicted output and measurement output. In addition, the research shown the formula to estimate the initial error covariance matrix P 0 based on using the inverse of the information matrix. Due to the use of batch data of N measurements, this method has some delay in the response to the change of sensor noise statistics. The delay time depends on the selection of the size of the estimation window N. Moreover, it also requires considerable computation and memory resources for the storage and processing of these N measurement times.
Macías and Expósito introduced a real-time adaptive method for self-tuning of the process noise covariance under the sudden changes of the input signal. 8) The Q matrix is updated in each time interval step using the norm of difference between the model predicted value and filter estimated value until the stopping condition is satisfied. However, the drawback of this research is that the Q matrix component may consist of measurement errors and the filter might be unstable if the measurement noise is very large compared with process noise.
The purpose of this paper is to propose an automatic and real-time tuning procedure for the Q and R matrices of the UKF in the redundant filter. This proposed algorithm is applied for a microsatellite named TSUBAME which is being developed at the Laboratory of Space Systems at the Tokyo Institute of Technology. In TSUBAME, various types of attitude sensors such as three-axis-magnetometer (TAM), Sun sensors, MEMS gyroscopes, star tracker (STT) system and fiber optic gyroscopes (FOG) are available. The main estimator uses a STT system and FOG while the redundant estimator only uses TAM and MEMS gyroscopes or FOG. The proposed filter uses the attitude data from the STT system in the main estimator as the truth reference data to compare with the estimated attitude of the redundant estimator. Based on the difference between the two attitudes above, the cost function is generated in each filter step. Then, it is minimized by a numerical optimization process using the downhill simplex algorithm to find the optimal values of the Q and R matrices. Through the tuning process, the convergence speed and estimated accuracy of the redundant estimator are improved. Based on that, to reduce the system power consumption, the proposed filter is used with the short duration and frequent repeat of turnon time of the main estimator data.
The organization of the paper is as follows. In section 2, the spacecraft rotation kinematics using quaternion representation and sensor models are briefly reviewed. In section 3, the UKF algorithm and the application of the UKF to satellite attitude estimation using attitude local error representation are presented. In section 4, the details of the real-time tuning unscented Kalman filter (RTUKF) algorithm are presented and discussed. In section 5, some simulation results are given with the discussions of stability, robustness, convergence speed and system power consumption. Finally, in section 6, the conclusions are given.

Attitude Kinematics and Sensor Models
In this section, a brief review of the attitude kinematics equation of motion using quaternions is shown. Then, the models of gyro and attitude sensor are briefly reviewed.

Attitude kinematics
The quaternion is defined by Ã T is the vector part and q 4 is the rotation part.
The quaternion representation is desirable because of its singularity free property. However, the norm constraint must be maintained. Since a four-dimension vector is used to represent three dimensions, the quaternion has a single constraint given by q T q ¼ 1. The attitude matrix is calculated as a quadratic function of q, that is where I 3Â3 is the 3 Â 3 identity matrix and ½&Â is the cross matrix defined as The quaternion kinematics differential equation is defined below. 6) Nov. 2013 H. X. LE and S. MATUNAGA: Real-Time Tuning Unscented Filter for a Redundant Attitude Estimator where ! is a 3 Â 1 vector of satellite angular velocity with respect to the initial frame and ÄðqÞ Equation (3) can be rewritten in the discrete form as and The subscript k denotes a measurement that is taken at time t k .

Sensor models
The satellite angular rate is measured by the gyro and its widely used model is given by 6) where! !ðtÞ is the measured angular velocity, !ðtÞ is the true angular velocity, ðtÞ is the gyro bias, and v ðtÞ and u ðtÞ are independent zero-mean Gaussian white-noise processes with standard deviations ' u and ' v . Gyro bias ðtÞ ¼ 1 2 3 Â Ã T is the unknown vector which is considered constant but can vary at each sampling step. The measurement vector sensor can be modeled as 6) y where v k $ N r ð0; R k Þ is a zero-mean Gaussian white noise with covariance R k , m is the number of vector measurements, b ik is a vector measurement in spacecraft body frame of ith sensor, q k is an attitude quaternion, Aðq k Þ is a corresponding attitude direction cosine matrix, r ik corresponds with the ith reference vector model in Earth-fixed frame and all values are represented at time t k .

UKF for Attitude and Bias Estimation
In this section the UKF is reviewed and then a UKF using a generalization of the Rodrigues parameters is derived for attitude and bias estimation.

Overview of the UKF
The filter is derived for a discrete-time nonlinear system modelled as where x k is an n-dimension state vector, u k is an input control, w k $ N r ð0; Q k Þ is a state process noise,ỹ y k is an observation vector of dimension m and v k $ N r ð0; R k Þ is a measurement noise vector. The process noise and measurement noise are assumed to be zero mean Gaussian noises with known covariance Q k and R k , respectively. The system function f and measurement function h with input timevalue-k can be time-varying. The sampling period is denoted by The idea of the UKF algorithm is to produce a set of well-chosen sampling points (sigma points) around the current estimated state based on its covariance. Then, these points are propagated through a nonlinear map to get an accurate estimation of the mean and covariance of the mapping results. Assume that the state vector x has the current estimated meanx x þ k and covariance P þ k . The set of 2n þ 1 sigma points, formed by the matrix 1 þ k , is generated by is the jth column of the square root matrix, ¼ ffiffiffiffiffiffiffiffiffiffiffi ffi n þ ! p , and ! ¼ 2 ðn þ Þ À n is the scaling parameter. The constant determines the spread of the sigma points aroundx x þ k and is usually set to a small positive value (e.g., 10 À4 10 À1 ). The constant is a secondary scaling parameter which is usually set to zero or 3 À n. The transformed set of sigma points through the nonlinear system is carried out for each point by The posterior mean and covariance are given bŷ The weights, W m i and W c i , are calculated using is used to incorporate prior knowledge of the distribution of state vector x (for Gaussian distributions, ¼ 2 is optimal 9) ).
The mean measurement and output covariance are given by The innovation covariance is calculated using Then, the cross-correlation matrix is determined by The Kalman gain matrix is approximated from the innovation covariance and cross-correlation matrix using Finally, updated state and covariance are calculated bŷ 3.2. Attitude and bias estimation using the UKF If the quaternion is directly used in the state vector of filters, the update process with additive approach may break the norm constraint condition. To solve this issue, we can use the local-error representation given by where q is the real attitude andq q is the estimated attitude.
The local-error used in the state vector is represented by a vector of generalized Rodrigues parameters as where a is a parameter from zero to one and f is a scale factor. Note that when a ¼ 0 and f ¼ 1, Eq. (30) gives the Gibbs vector; when a ¼ f ¼ 1, Eq. (30) gives the standard vector of modified Rodrigues parameters. The inverse transformation from p to q is given by Since the three-component of generalized Rodrigues parameters, in which the singularity can be placed anywhere from AE180 to AE360 6,10) is used to represent small attitude errors, the singularity hardly occurs in practice.
The state vector of UKF is chosen aŝ where p p þ k is the local error of spacecraft attitude represented in generalized Rodrigues parameters, þ k is the gyro bias vector and n is the dimension of state vector.
(2) Time update: Each sigma point is transformed sequentially through the system function f and the measurement function h to obtain the transformed samples. The set of current estimated quaternion is calculated by where q þ k ðiÞ is converted from 1 p þ k ðiÞ. The estimated angular velocities are calculated from the measured angular ratẽ ! ! k by!
The quaternions of local errors are calculated from the updated quaternion q À kþ1 ðiÞ by All sigma points are now updated by converting q À kþ1 ðiÞ to 1 p À kþ1 ðiÞ and setting 1 À kþ1 ðiÞ ¼ 1 þ k ðiÞ because gyros bias is assumed to be the constant number. Then, the posterior mean and covariance can be calculated by Eqs. (16) and (17), respectively.
(3) Measurement update: Similarly, the mean measurement and output covariance are calculated by Eqs. (21) and (22) with Then, the innovation covariance and cross-correlation matrix are calculated using Eqs. (24) and (25), respectively. The Kalman gain matrix is computed by Eq. (26).
Equations (27) and (28) are used to update the state vector x x þ kþ1 and the error covariance matrix P þ kþ1 . Note in Eq. (27) thatỹ y kþ1 is the measurement value of the attitude sensor at the current time. After convertingx x p þ kþ1 to q þ kþ1 , the quaternion is updated bŷ Finally, the local error is reset to zero for the next step, sô The state vectorx x þ kþ1 , attitude quaternionq q þ kþ1 and error covariance matrix P þ kþ1 are updated and ready for the next loop.

Real-Time Tuning UKF for a Redundant Attitude Estimator
In this section, the details of the proposed method, named the RTUKF algorithm, are presented and discussed. As mentioned above, the proposed method is applied for the redundant estimator, which supports the working of the main estimator. Moreover, in the case the satellites suddenly enter the power saving mode, if the redundant estimator Nov. 2013 H. X. LE and S. MATUNAGA: Real-Time Tuning Unscented Filter for a Redundant Attitude Estimator starts at this time, it will take time to converge to the steady estimated value. Therefore, it is preferred that the redundant estimator, which uses low-power, low-accuracy sensors such as TAM, remains on. The main estimator has a higher level of estimated accuracy and stability because of the use of high-power, high-accuracy sensors such as STT and FOG. In order to improve the quality of this redundant estimator, the estimated results from the main estimator are used as the truth reference data for the tuning process. The tuning process finds the optimal values of the Q and R matrices for each calculation step by minimizing the cost function. The cost function is calculated based on the truth reference value from the main estimator so the tuning process can only be implemented when the estimated result of the main estimator exists. The general algorithm of the RTUKF is shown in Fig. 2. Depending on the availability of reference data, the current optimal values Q Ã and R Ã or the constant optimal values Q 0 and R 0 are selected for the conventional UKF. Since RTUKF is implemented in each calculation step as a real-time algorithm, the tuning parameters are strongly affected by the noise of sensors which are used for the redundant estimator.

Parameterization and optimal function
It takes a large amount of calculation effort if all parameters of the Q and R matrices are tuned. In conventional UKF, 5,6,11) Q and R are defined as the constant matrices using static optimal values as where ' TAM is the standard deviation of TAM noise and ' u and ' v are introduced in the gyro model part.
In this proposed method, Q and R are tuned by two positive scalar scaling factors Q and R as Q ¼ Q Q 0 and R ¼ R R 0 . To ignore the positive number constraint, Q , R can be calculated through $ Q and $ R by where $ Q , $ R are arbitrary numbers in R. The cost function is defined as where P pp kþ1 2 R 3Â3 is a sub-matrix of the error covariance part matrix P kþ1 2 R 6Â6 at the step k defined as and p kþ1 is a generalized Rodrigues parameters representation of the propagated error quaternion, which is calculated as q k ¼ q STT q À1 est at the step k. q STT is the attitude quaternion output from the STT system and q est is the estimated attitude quaternion from the UKF output in the redundant estimator with the current values of the Q and R matrices.

Numerical optimal algorithm
The real-time tuning process should include a numerical optimization algorithm to find the best match of the Q and R matrices in real-time for each attitude estimation step. The downhill simplex method 12) is chosen as the optimal algorithm for the tuning process. This is a multidimensional optimization method which uses geometric relationships to aid in finding function minimums. One distinct advantage of this method is that it does not require the derivative of the function. Therefore it works well with a strongly nonlinear cost function. The downhill simplex creates its own pseudoderivative by evaluating enough points of the function to define the derivative for each independent variable. This method revolves around a simplex: a geometric object that contains n þ 1 vertices where n is the number of independent variables. For instance, in the one-dimensional optimization, the simplex contains two points and represents a line segment. In the two-dimensional and threedimensional optimizations, the simplex is a three pointed triangle and a tetrahedron, respectively. The function to be minimized is evaluated at each of the N þ 1 simplex vertices.
As the RTUKF needs to tune two parameters, the downhill simplex has three vertexes. Each vertex has a function value that is exactly same as the cost function defined in Eq. (43). The goal of this algorithm is then to move the simplex away from points with larger value. The downhill simplex algorithm can employ several methods for moving the simplex downhill. These methods include reflection, expansion, contraction and reduction. The chosen one depends on both data known from the current simplex location as well as data from previous moves.
To implement the downhill simplex for RTUKF, it is necessary to convert all equations which include UKF formulas and the calculation of the cost function based on estimated error p kþ1 into the equivalence function. This equivalence function also returns the same scalar value of the cost function as in Eq. (43). Figure 3 shows the downhill simplex algorithm inside the ''optimal loop of tuning process'' block in Fig. 2. This loop is repeated until the difference between the maximum and minimum of cost functions at all vertexes smaller than a predefined number ".

Numerical Simulation Results
In this session, the proposed method is applied for a microsatellite named TSUBAME. We consider three types of attitude sensors of the TSUBAME including TAM, STT and FOG. TAM remains on in all working modes and it is used as the main vector sensor in the redundant estimator. STT is assumed to directly provide the estimated attitude of the satellite in the quaternion representation with the accuracy less than 30 arcseconds right after turned on. Main specifications of all sensors are shown in Table 1.
The simulation process is setup as in Fig. 4. The satellite orbit is generated from the TLE file based on the SGP4 model, and then the reference Earth magnetic field in the ECI frame is calculated using the IGRF-11 model. From satellite dynamics and kinematics models, the satellite attitude and satellite angular velocity are generated. Through the gyro model and TAM model, the omega measured on FOG! ! k and TAM dataỹ y kþ1 are created as the simulated sensor data including several types of noise. STT data is also generated by adding the zero-mean white Gaussian noise to the real quaternion.
All simulation scenarios are implemented as the Monte Carlo simulation with 50 runs in which the initial attitude error and gyro bias on each rotation axis are randomly generated alternately within AE30 deg and AE90 (deg/h). The initial error covariance matrix is P 0 1e À3 I 6Â6 and the initial state vector is x 0 0 ... 0 ½ T . Total simulation time is 2 h with sampling frequencies of all sensors set at 1 Hz. The plotted results regarding the attitude and the bias estimated errors are the average values from all 50 runs. Figure 5 shows the first simulation scenario results where the first four curves show the norm of the RTUKF estimated errors corresponding to four different cost functions defined as The fifth curve is the conventional UKF using additional data from STT as the measurement vector-sensor. The sixth curve is the conventional UKF using only TAM and FOG data. The seventh curve plotted in the hidden line shows the availability of STT. Note that only when STT exists, the tuning process of RTUKF can be implemented. From simulation results, the first curve is the best one in terms of convergence speed and stability. The first candidate converges to STT error level within 100 s while the second and third converge within 200 and 250 s. The fourth curve shows that if there is no extra reference data from STT, the optimal process does not produce any improvement. It is also proved that KFs themselves are already the optimal solutions. Figure 6 shows details of the first three values in the diagonal vector, diag ðPÞ ' 2 11 . . . ' 2 66 Â Ã T of the error covariance matrix P corresponding to the first and third candidate cost functions. The first cost function has a big step of con-vergence when RTUKF is applied. This is one of the causes explaining the fastest convergence speed of the estimated error curve shown in Fig. 5. Figure 7 presents the value of the tuning parameter $ R corresponding to the first, third and fourth candidate cost functions in the same simulation. As mentioned above, the tuning parameter is very frequently changed to resist the effect of the TAM noise.
The second simulation scenario, for which results are shown in Fig. 8, is designed to verify the effect of the types of noise in TAM. There are three types of noise distribution (normal, modified-zero-mean Poisson and uniform distribution) that are taken into account. In this figure, the effect of noise distribution type is not clear during the normal working condition (without the reference data from STT). The normal distribution is only slightly better compared with others during the recovery time to the normal condition after the tuning process of RTUKF is turned off.
The power generation of microsatellites is very limited so the satellites may sometimes have the critical power consumption condition. In this case, they need to reduce their consumption power as soon as possible. However, the mission of the satellites still requires a good enough estimated Time (s) 0.5*dp'*inv(P(1:3,1:3))*dp 0.5*dp'*dp 0.5*dp'*(P(1:3,1:3))*dp 0.5*(y-y1)'*inv(P(1:3,1:3))*(y-y1) UKF using STT data Original UKF STT availability  (0.5*dp'*inv(P(1:3,1:3))*dp) (0.5*dp'*(P(1:3,1:3))*dp) (0.5*(y-y1)'*inv(P(1:3,1:3))*(y-y1)) Parameter of R matrix for RTUKF  accuracy which the redundant system cannot only archive. To solve this problem, the RTUKF is used to reduce the system power consumption but still keep the estimated accuracy at a good enough level. At that time, the RTUKF is used with the short duration and frequent repeat of turnon time of the main estimator data, the STT data. In the last simulation scenario, a series of simulation samples is designed with a constant duration working time of STT (100 s) in the variable periods of time, such as 1,000, 667, 500, 250 and 200 s corresponding to 10, 15, 20, 40 and 50% using time of STT, respectively. In these simulation samples, the UKF employing STT as the measurement vector-sensor are compared with the RTUKF. Note that to calculate the power consumption, the STT system has two separate STT units with 1.6 W power consumption on each one, while the power consumption of TAM is 0.8 W for all three axis measurements. The chart in Fig. 9 shows the summary of estimated errors and the power consumption in different working conditions. The estimated errors are computed at the highest point in the last turn on time of STT (see details in Fig. 10). Considering the cases of 10 and 50% using time of STT, respectively, the attitude error of ''UKF using STT data'' method is reduced 43.8% while the attitude error of RTUKF method is reduced 80.8%. It shows the strong advantage of the RTUKF in terms of convergence speed compared with the convenient UKF when the high-accuracy reference data only exists for a short time. From this characteristic, depending on the percentage of usage time of STT, the estimated error of the RTUKF can be archived from 0.125 deg without STT reference data up to 0.0083 deg with full-time-usage of STT.
In more detail, regarding the simulation sample in the case of 20% usage time of STT, the attitude estimated error and the gyros bias estimated error are plotted in Figs. 10 and 11 respectively. The RTUKF only has the reference data on the satellite attitude and the cost function to consider the attitude error. Therefore, it affects the estimated gyros bias values, and their estimated error are slightly increased during the first optimal loops of the tuning process. However, this error is still acceptable.

Conclusion
This paper presented a new filter algorithm named the RTUKF for a redundant system. The RTUKF uses a downhill simplex search to tune the process noise and measurement noise covariance matrices of the conventional UKF in each time step. The proposed filter parameters were tuned in real-time by minimization of the cost function calculated from the difference between the filter estimated attitude in the redundant estimator and the reference attitude in the main estimator. The RTUKF had higher accuracy and faster convergence speed than the conventional UKF employing STT as the measurement vector-sensor. With the advantage on convergence speed, the proposed estimator can be used with the suitable duration and repeat frequency of turn-on time of the main estimator to meet mission requirements on either the estimated accuracy or the system power consumption. Illustrated by the numerical simulation results, it is evident that there always exist the optimal values of tun-