Simplex Back Propagation Estimation Method for Out-of-Sequence Attitude Sensor Measurements

For a small satellite, the processor onboard the attitude determination and control system (ADCS) is required to monitor, communicate, and control all the sensors and actuators. In addition, the processor is required to consistently communicate with the satellite bus. Consequently, the processor is unable to ensure all the sensors and actuators will immediately respond to the data acquisition request, which leads to asynchronous data problems. The extended Kalman filter (EKF) is commonly used in the attitude determination process, but it assumes fully synchronous data. The asynchronous data problem would greatly degrade the attitude determination accuracy by EKF. To minimize the attitude estimation accuracy loss due to asynchronous data while ensuring a reasonable computational complexity for small satellite applications, this paper proposes the simplex-back-propagation Kalman filter (SBPKF). The proposed SBPKF incorporates the time delay, gyro instability, and navigation error into both the measurement and covariance estimation during the Kalman update process. The performance of SBPKF has been compared with EKF, modified adaptive EKF (MAEKF), and moving–covariance Kalman filter (MC-KF). Simulation results show that the attitude estimation error of SBPKF is at least 30% better than EKF and MC-KF. In addition, the SBPKF’s computational complexity is 17% lower than MAEKF and 29% lower than MC-KF.


Introduction
The development and manufacturing of the small satellite have been a focus in the industry over the past decades. It is reported that approximately 350 small satellites with a mass of less than 200 kg have been launched in 2021 [1]. The Satellite Technology And Research Centre (STAR) at the National University of Singapore is currently developing three satellites named Lumelite-1 to -3 for formation flying programs [2]. Each satellite has a wet mass of 18 kg. The Lumelite satellites are planned to be launched in 2023.
The Lumelite satellite bus system primarily consists of the attitude determination and control subsystem (ADCS), the onboard computer (OBC) system, the electrical power system (EPS), and the communication interface module (CIM). Each of these subsystems is required to continuously monitor its electrical status, communicate with the sensor or actuator that connects to the subsystem, and communicate with OBC. For example, the ADCS's digital signal processor (DSP) is required to ensure a stable power supply to all sensors, and have no over-current drawn by any actuator, while continuously requesting data from the sensor and sending a command to the actuator at a fixed sampling interval. Due to the reason that each sensor requires a certain processing time upon receiving the request, and the DSP is processing other tasks concurrently the updated data from each sensor and actuator would be received at different time instances, which results in unsynchronized-sensor and actuator information.
Estimation algorithms, such as the extended Kalman filter (EKF) [3], Unscented Kalman filter (UKF) [4,5], Cubature Kalman filter (CKF) [6] and particle filter [7], are Figure 1 illustrates the input and output (I/O) of ADCS between sensors, actuators, satellite bus systems, memory, payload, and DSP. Three Universal Asynchronous Receiver/Transmitter (UART) interfaces are used for the communication between DSP and magnetometer, communication between DSP and GPS receiver, and ADCS debugging purposes. The GPS receiver also directly provides a 1 pulse-per-second (1PPS) signal to DSP via a dedicated interface. The DSP communicates with the satellite bus system via a Controller Area Network (CAN) interface. The additional CAN interface is used for propulsion system communication purposes. The DSP controls the magnetic torquer via pulse-width modulation (PWM) interface, while it receives the coarse sun sensor data via an analog-to-digital converter (ADC) I/O. Lastly, the DSP communicates with a fine sun sensor and reaction wheel via RS-485 and RS-232 interfaces, respectively.

Tasks of Attitude Determination and Control System
algorithm. The proposed SBPKF is presented in Section 5. Section 6 presents the s tion and results. Finally, Section 7 concludes the paper. Figure 1 illustrates the input and output (I/O) of ADCS between sensors, act satellite bus systems, memory, payload, and DSP. Three Universal Asynchrono ceiver/Transmitter (UART) interfaces are used for the communication between D magnetometer, communication between DSP and GPS receiver, and ADCS deb purposes. The GPS receiver also directly provides a 1 pulse-per-second (1PPS) si DSP via a dedicated interface. The DSP communicates with the satellite bus syste Controller Area Network (CAN) interface. The additional CAN interface is used f pulsion system communication purposes. The DSP controls the magnetic torq pulse-width modulation (PWM) interface, while it receives the coarse sun sensor d an analog-to-digital converter (ADC) I/O. Lastly, the DSP communicates with a f sensor and reaction wheel via RS-485 and RS-232 interfaces, respectively. The I/O diagram in Figure 1 indicates that the DSP is required to perform v tasks to ensure the stability of the satellite. Furthermore, Table 1 indicates the sa period of each task that must be handled by the DSP. During the early phase of fir development and testing, it had been noticed that the sensor data acquisition pro quires at least 55 milliseconds to 75 milliseconds (represented by , and , in 2) per cycle. As shown in Figure 2, both the sun sensor and magnetometer require a time to respond to the request sent by DSP ( , and , , respectively) when the task is initiated. For the ideal scenario, , and , shall be less than 5 millisecon noted that the test process only involves both ADC tasks and sensor tasks. It is ex that the sensor data acquisition process requires much longer time in the full fir implementation.  The I/O diagram in Figure 1 indicates that the DSP is required to perform various tasks to ensure the stability of the satellite. Furthermore, Table 1 indicates the sampling period of each task that must be handled by the DSP. During the early phase of firmware development and testing, it had been noticed that the sensor data acquisition process requires at least 55 milliseconds to 75 milliseconds (represented by t S,r and t B,r in Figure 2) per cycle. As shown in Figure 2, both the sun sensor and magnetometer require a certain time to respond to the request sent by DSP (t S,r and t B,r , respectively) when the sensor task is initiated. For the ideal scenario, t S,r and t B,r shall be less than 5 milliseconds. It is noted that the test process only involves both ADC tasks and sensor tasks. It is expected that the sensor data acquisition process requires much longer time in the full firmware implementation.  In Figure 2, the attitude determination (AD) task and sensor task are performed at different frequencies. In addition to the delay in the sensor data acquisition process, the asynchronization issue between sensor data and the AD task (or by sun sensor and by magnetometer, respectively) would degrade the overall attitude determination accuracy. Thus, an algorithm to minimize the accuracy degradation due to the sensor data synchronization issue is required. As previously mentioned, the sun sensor or magnetometer measurement will only be considered as a delayed measurement if either or is more than five milliseconds. The five milliseconds is the expected response time of the sensor under the ideal operation scenario.

Extended Kalman Filter
The small satellite is typically equipped with sun sensors, magnetic torquers, and microelectromechanical systems (MEMS) based inertial measurement unit (IMU). The MEMS IMU device includes both gyroscope and magnetometer sensors. The output of the MEMS IMU gyroscope is given by [24] ̃k = + + g + β (1) where overhead notation of "." denotes measurement, ̃k denotes the measured threeaxis body rate, denotes the truth three-axis body rate, g denotes the output noise of the gyroscope, denotes the truth gyroscope bias, and β denotes the gyroscope bias due to the angular random walk and bias instability.
The general EKF process consists of gain matrix computation, estimated state and covariance update, and estimated state and covariance propagation. For EKF-based attitude determination, the objective is to estimate both the satellite's attitude (or quaternion) and gyroscope bias. Thus, the estimated state vector is ̂k = [̂k̂T]. The k is the vector component of the quaternion vector, , such that [25] where the overhead notation of "." denotes the estimated value, and ̂4 , is the scalar of the quaternion vector with ̂4 , = √1 −̂k T̂k .
In Figure 2, the attitude determination (AD) task and sensor task are performed at different frequencies. In addition to the delay in the sensor data acquisition process, the asynchronization issue between sensor data and the AD task (or δt S by sun sensor and δt B by magnetometer, respectively) would degrade the overall attitude determination accuracy. Thus, an algorithm to minimize the accuracy degradation due to the sensor data synchronization issue is required. As previously mentioned, the sun sensor or magnetometer measurement will only be considered as a delayed measurement if either δt S or δt B . is more than five milliseconds. The five milliseconds is the expected response time of the sensor under the ideal operation scenario.

Extended Kalman Filter
The small satellite is typically equipped with sun sensors, magnetic torquers, and microelectromechanical systems (MEMS) based inertial measurement unit (IMU). The MEMS IMU device includes both gyroscope and magnetometer sensors. The output of the MEMS IMU gyroscope is given by [24] ω k = ω k + β + η g + η β (1) where overhead notation of " ." denotes measurement, ω k denotes the measured three-axis body rate, ω k denotes the truth three-axis body rate, η g denotes the output noise of the gyroscope, β denotes the truth gyroscope bias, and η β denotes the gyroscope bias due to the angular random walk and bias instability. The general EKF process consists of gain matrix computation, estimated state and covariance update, and estimated state and covariance propagation. For EKF-based attitude determination, the objective is to estimate both the satellite's attitude (or quaternion) and gyroscope bias. Thus, the estimated state vector isx k = ρ T k β T . The ρ k is the vector component of the quaternion vector, q, such that [25] where the overhead notation of "." denotes the estimated vue, andq 4,k is the scalar of the quaternion vector withq 4,k = 1 − ρ T k ρ k . Given a pair of measurement vectors from sensors, and their corresponding estimated measurement vector, bothq k and β are updated via the following process [25] where Ξ(q k ) is defined in [25] and In (4), y ≡ y T S y T B T , where the subscript "S" denotes the sun vector and subscript "B" denotes the earth's magnetic field vector. In addition, the superscript "−" denotes pre-update, superscript "+" denotes post update, and K denotes the Kalman filter gain matrix [26] where R k denotes measurement noise covariance and H k denotes the Jacobian matrix of measurement vectors. Moreover, the mth sensor component of the estimated measurement vector in (4) can be expressed as [25] where A q − k denotes the estimated attitude matrix at time t k , andê m (a unit vector) is the estimated mth sensor's measurement in the earth center inertia (ECI) reference frame.
It is noted that the derivation ofê S associated toŷ S is detailed in [27], andê B is provided in Appendix A and (31). In addition, H k in (5) is given as [26] For standard EKF, ∂ŷ S /∂ρ, ∂ŷ B /∂ρ, ∂ŷ S /∂β and ∂ŷ B /∂β are defined as [28] ∂ŷ Here, [a×] is the cross-product matrix [25]. By assuming that there is no crosscorrelation between the sun vector and the earth magnetic field vector, R k can be written as where R S,k is the covariance associated with the sun sensor, and R B,k is the covariance associated with the magnetometer. In general, the measurement vector, b m output by the m th sensor in the satellite body reference frame always contains an error, with the standard deviation of ν m : where b m denotes the truth measurement reading of the m th sensor. Furthermore, b m maybe not necessarily a unit vector (e.g., Earth's magnetic field vector). Therefore, y m represents the normalized vector of b m such that y m = b m / b m , where . denotes the vector's magnitude. Then, the general expression of R m,k in (9) with the consideration of vector normalization is given as [29]: with E{.} denotes expected value, trace R m,k denotes the summation of diagonal elements of R m,k . In addition, the state error covariance update is given as [26] Using the following definition,ω k = ω k − β, the nonlinear quaternion propagation model is given as [30] .
where Ω ω k is defined in [30]. Due to the nature of quaternion multiplication, the integration of (14) is often highly complex. Instead, the quaternion propagation can be simplified using the N x -step summation model [31]: The state error covariance propagation model is given by [26] .
where [25] The standard EKF shows that the state vector update process in (3) to (5), and the Jacobian matrix in (8) do not consider the unsynchronized sensor data problem. On the other hand, the proposed SBPKF and MC-KF include the sensor data delay and its associated sensitivity matrix and error covariance in the measurement model estimation. The MC-KF algorithm is presented in the next Section.

Moving Covariance Kalman Filter
The presented MC-KF in this section is based on [23,32]. It has a similar estimation process as DV-DTCF in [12] but a different state vector fusion approach. Let us define the following time instance, t m = t k − δt m , where subscript "m" is a general representation of either sun sensor,"S" or magnetometer "B". First, given the sun sensor and magnetometer data arrive at different time instances, there are two sets of state correction vector, ∆x m,t m where K m is the Kalman filter gain matrix that is derived with respect to the sensitivity matrix of each sensor: where R m,k is defined in (11).
From (18), two sets of updated state vectors,x + m, t m (orx + S, t S andx + B, t B ) and two sets of state error covariances, P + m, t m (or P + S, t S and P + B, t B ) will be obtained by using a similar formulation in (3) and (13). Next, bothx + S, t S , P + S, t S andx + B, t B , P + B, t B pairs are required to be synchronized and smoothed. First, let's define the following time instance, t u based on the following sensor delay information As presented in Figure 3, the delay of 5 ms is assumed to be the minimum time required by a sensor to instantly reply to the sensor data request by ADCS DSP. Then, thê x P, t p and P P, t p pair that has a lesser delay with respect to t k , are propagated to t û  For generalization purposes, the sensor data with a higher delay is labeled as w th sensor. Based on [23], we define the following matrices, P/W and P−W For generalization purposes, the sensor data with a higher delay is labeled as w th sensor. Based on [23], we define the following matrices, P P/W and P P−W where K m,t u and H m,t u are Kalman filter gain and sensitivity (derivation see (7), (8), (33) and (34)) matrices at t u . To ensure the positive definite of P P−W matrix, we have In (25), "diag(A)" denotes the matrix only contains the diagonal elements of matrix A. The purpose is to ensure the matrix P * P−W is invertible. Then, the smoothed state vector is given asx The quaternion update corresponding to (26) is given in (3). On the other hand, the quaternion subtraction to compute δx t u is given as follows: with Ψ(q) is defined in [25]. Finally, the covariance is smoothed as follows Lastly, both state vectors,x + t u and covariance, P + t u , are propagated to the next time step, x − k+1 and P − k+1 via (15) and (16).

Simplex-Back-Propagation Kalman Filter
The proposed SBPKF uses the similar Kalman filter process as EKF, but with different measurement and covariance models. The SBPKF has taken the delay between sensor data acquisition time and Kalman sampling time, δt m (or known as δt S and δt B in Figure 2). As such, the delayed measurement vector for both the sun sensor and magnetometer in (6) is rewritten as: It is noted that the Earth's magnetic field vector in the ECI reference frame, e B is computed based on the satellite position and velocity at t k − δt B where A ECI NED denotes the transformation matrix of NED to ECI reference frame, and A ECI NED will not be shown in this paper, but the associated algorithm and Matlab code are available in [33]. In addition, r t k −δt B and v t k −δt B can be approximated as From (30), the measurement vector is in terms of both quaternion and gyro bias. Therefore, both ∂ŷ S /∂β and ∂ŷ B /∂β are no longer zero matrices. Instead, with the additional summation terms, (8) becomes The general expression of R m,k remains the same as in (11). However, (31) indicates the satellite positioning error contributes to the Earth's magnetic field modeling error. In addition, both η g and η β contribute additional errors in (30). Therefore, both R S,k and R B,k , which were originally formulated in (12) have been modified to become: It is noted that the term σ B r t k −δt B denotes additional error covariance due to the satellite position error: In (38), ∂B/∂h LLA denotes the partial derivative of B r t k −δt B , v t k −δt B in (31) with respect to longitude, φ, latitude, θ and radius, r (or h LLA = θ φ r T ). In addition, ∂h LLA /∂r ECEF denotes the partial derivative of conversion from position in the Earthcentered-Earth-fix (ECEF) reference frame to longitude, latitude, and radius. The formulation for ∂B/∂h LLA is available in Appendix A, and the formulation for conversion of ECEF position to longitude, latitude, and radius is detailed in [34].
The overall process flow for SBPKF (and MC-KF/MAEKF) is presented in Figure 3. When the satellite exits from the eclipse, the first pair of sun and Earth magnetic field vectors will be input into the quaternion estimation (QUEST) algorithm [35] to provide an initial quaternion vector, with the assumption of fully synchronized measurement vectors. The quaternion output from QUEST also guarantees the stability of EKF and SBPKF algorithms as the EKF-based algorithm is often susceptible to initial condition error [36]. As discussed in Section II, the respective measurement will only be considered as a delayed measurement for the filtering process if δt m > 5 ms. Thus, if both δt S and δt B are less than a given threshold (i.e 5 ms), standard EKF will be conducted. Otherwise, the delayed measurement vector in (30) will be computed. Once the measurement covariances are computed using (11), (35) to (37), the SBPKF updates the estimated states and state error covariance using the (3) and (13).

Simulation and Results
Monte Carlo simulations have been conducted to compare the quaternion and gyro bias estimation accuracy for the proposed SBPKF, EKF, MAEKF in [22], and MC-KF, with respect to sensor delay. The Monte Carlo simulation environment is illustrated in Figure 4. Three ADCS tasks are considered in the simulation, which are the sensor task, AD task, and AC task. The sensor task is simulated at a sampling rate of 200 ms, with selected sensor delays. The sensor delay to be range from 65 ms to 145 ms, with its standard deviation given in Table 2. The AC task simulates the truth quaternion vector that is corresponding to the satellite pointing profile in Figure 5. The AD task comprises the attitude determination algorithm, such as SBPKF, EKF, MAEKF, and MC-KF. During the eclipse period (or without sunlight), the satellite attitude control enters the momentum hold condition. The default attitude control of the satellite during the sunlight condition is the sun-pointing mode. The satellite performs nadir pointing when the angle between the nadir axis and sun vector in the satellite body frame is within the sun sensor's half-cone field-of-view (or 60 degrees). In Figure 5, the satellite begins to perform nadir pointing 15 min after entering the sunlight condition, for approximately 35 min.
vectors will be input into the quaternion estimation (QUEST) algorithm [35] to provide an initial quaternion vector, with the assumption of fully synchronized measurement vectors. The quaternion output from QUEST also guarantees the stability of EKF and SBPKF algorithms as the EKF-based algorithm is often susceptible to initial condition error [36]. As discussed in Section II, the respective measurement will only be considered as a delayed measurement for the filtering process if > 5 ms. Thus, if both S and B are less than a given threshold (i.e 5 ms), standard EKF will be conducted. Otherwise, the delayed measurement vector in (30) will be computed. Once the measurement covariances are computed using (11), (35) to (37), the SBPKF updates the estimated states and state error covariance using the (3) and (13).

Simulation and Results
Monte Carlo simulations have been conducted to compare the quaternion and gyro bias estimation accuracy for the proposed SBPKF, EKF, MAEKF in [22], and MC-KF, with respect to sensor delay. The Monte Carlo simulation environment is illustrated in Figure  4. Three ADCS tasks are considered in the simulation, which are the sensor task, AD task, and AC task. The sensor task is simulated at a sampling rate of 200 ms, with selected sensor delays. The sensor delay to be range from 65 ms to 145 ms, with its standard deviation given in Table 2. The AC task simulates the truth quaternion vector that is corresponding to the satellite pointing profile in Figure 5. The AD task comprises the attitude determination algorithm, such as SBPKF, EKF, MAEKF, and MC-KF. During the eclipse period (or without sunlight), the satellite attitude control enters the momentum hold condition. The default attitude control of the satellite during the sunlight condition is the sunpointing mode. The satellite performs nadir pointing when the angle between the nadir axis and sun vector in the satellite body frame is within the sun sensor's half-cone fieldof-view (or 60 degrees). In Figure 5, the satellite begins to perform nadir pointing 15 min after entering the sunlight condition, for approximately 35 min.      Table 2. It is noted that the right ascension of ascending node, the argument of perigee, and the initial true anomaly are randomly generated in each Monte Carlo simulation.
The sensors and gyro noises, GPS accuracy, sensor, and EKF sampling rate are listed in Table 2. The sun sensor, GPS, gyro, and magnetometer noises are based on the commercial off-the-shelf product, with the consideration of noise density and sensitivity. In addition, the configuration of MAEKF is based on details provided in [22]. Figure 6 compares the quaternion and gyro bias estimation error between the SBPKF, EKF, MAEKF, and MC-KF with respect to the estimation results without the sensor delay scenario. The average quaternion error magnitude and gyro bias error magnitude for the zero sensor delay scenario is 0.145 degrees and 2.8 millidegrees, respectively. The quaternion error ratio, ̿, and gyro bias error ratio, ̿ , in Figure 6 are given by:

Accuracy Comparison
where M d is the total number of samples for the delayed sensor scenario, M 0 is the total number of samples for no sensor delay scenario, ∆ k is the quaternion error, and ∆ k is For each Monte Carlo simulation, the performance of SBPKF, EKF, MAEKF, and MC-KF are evaluated for one sunlight period (or approximately 60 min). The satellite orbital parameters are provided in Table 2. It is noted that the right ascension of ascending node, the argument of perigee, and the initial true anomaly are randomly generated in each Monte Carlo simulation.
The sensors and gyro noises, GPS accuracy, sensor, and EKF sampling rate are listed in Table 2. The sun sensor, GPS, gyro, and magnetometer noises are based on the commercial off-the-shelf product, with the consideration of noise density and sensitivity. In addition, the configuration of MAEKF is based on details provided in [22]. Figure 6 compares the quaternion and gyro bias estimation error between the SBPKF, EKF, MAEKF, and MC-KF with respect to the estimation results without the sensor delay scenario. The average quaternion error magnitude and gyro bias error magnitude for the zero sensor delay scenario is 0.145 degrees and 2.8 millidegrees, respectively. The quaternion error ratio, = q, and gyro bias error ratio, = β, in Figure 6 are given by:

Accuracy Comparison
where M d is the total number of samples for the delayed sensor scenario, M 0 is the total number of samples for no sensor delay scenario, ∆q k is the quaternion error, and ∆β k is gyro bias error at k th sample with subscript "d" representing delayed sensor, scenario, and the subscript "0" representing no sensor delay scenario, respectively. gyro bias error at k th sample with subscript "d" representing delayed sensor, scenario, and the subscript "0" representing no sensor delay scenario, respectively.
(a) (b)  Figure 6a shows that the SBPKF has the lowest quaternion error magnitude, followed by MAEKF. The SBPKF's quaternion estimation error ratio is approximately 189% (or 0.274 degrees) when compared to the no-sensor-delay scenario. The MAEKF quaternion error ratio is approximately 200% (or 0.29 degrees). Figure 6a also shows that without the simplex-back-propagation method, the quaternion estimation error ratio of EKF is increased to 220-250% (or 0.319 to 0.363 degrees), and the quaternion estimation error ratio of MC-KF is increased to 250-300% (or 0.363 to 0.435 degrees). Figure 6b shows that the EKF has the lowest gyro bias error magnitude. However, the gyro bias error magnitude difference between EKF and the SBPKF can be considered negligible as both are in the range of millidegrees (or approximately three millidegrees). On the other hand, the MAEKF gyro bias error is three times higher than the SBPKF, even though it has a similar quaternion estimation error as the SBPKF. The additional error that occurs in MAEKF could be due to the reason that the MAEKF assumes all sensor data received at the same time instance. However, in practice, each sensor has a different time delay.
Although the MC-KF is designed to compensate for the error due to data delay, the introduction of unknown gyro bias as the initial condition has greatly degraded the performance of the MC-KF. From Figure 6b, the results show that the MC-KF has a large gyro bias error ratio as compared to the SBPKF and EKF. The MC-KF's gyro bias error is approximately five times higher than the scenario without data delay. The MC-KF uses one measurement vector to update each ̂S ,t S and ̂B ,t B (see (18) and (26)). However, the attitude determination process requires at least a pair of measurement vectors to effectively estimate the quaternion vector. Furthermore, the additional covariance update process in (29) causes the state error covariance to be underestimated in MC-KF. The study in [37] has shown that an additional covariance update process within an iteration without a proper smoothing procedure could degrade the estimation accuracy. Therefore, the results in Figure 6 show that MC-KF would require a more complex derivation and implementation to ensure a stable estimation performance.
On the other hand, the proposed SBPKF considers the time delay of sensor data, δ m when estimating the measurement vector in (30). In addition, (35) to (38) consider the  Figure 6a shows that the SBPKF has the lowest quaternion error magnitude, followed by MAEKF. The SBPKF's quaternion estimation error ratio is approximately 189% (or 0.274 degrees) when compared to the no-sensor-delay scenario. The MAEKF quaternion error ratio is approximately 200% (or 0.29 degrees). Figure 6a also shows that without the simplex-back-propagation method, the quaternion estimation error ratio of EKF is increased to 220-250% (or 0.319 to 0.363 degrees), and the quaternion estimation error ratio of MC-KF is increased to 250-300% (or 0.363 to 0.435 degrees). Figure 6b shows that the EKF has the lowest gyro bias error magnitude. However, the gyro bias error magnitude difference between EKF and the SBPKF can be considered negligible as both are in the range of millidegrees (or approximately three millidegrees). On the other hand, the MAEKF gyro bias error is three times higher than the SBPKF, even though it has a similar quaternion estimation error as the SBPKF. The additional error that occurs in MAEKF could be due to the reason that the MAEKF assumes all sensor data received at the same time instance. However, in practice, each sensor has a different time delay.
Although the MC-KF is designed to compensate for the error due to data delay, the introduction of unknown gyro bias as the initial condition has greatly degraded the performance of the MC-KF. From Figure 6b, the results show that the MC-KF has a large gyro bias error ratio as compared to the SBPKF and EKF. The MC-KF's gyro bias error is approximately five times higher than the scenario without data delay. The MC-KF uses one measurement vector to update eachx S, t S andx B, t B (see (18) and (26)). However, the attitude determination process requires at least a pair of measurement vectors to effectively estimate the quaternion vector. Furthermore, the additional covariance update process in (29) causes the state error covariance to be underestimated in MC-KF. The study in [37] has shown that an additional covariance update process within an iteration without a proper smoothing procedure could degrade the estimation accuracy. Therefore, the results in Figure 6 show that MC-KF would require a more complex derivation and implementation to ensure a stable estimation performance.
On the other hand, the proposed SBPKF considers the time delay of sensor data, δt m when estimating the measurement vector in (30). In addition, (35) to (38) consider the additional propagation error due to the measurement model used in (30). Thus, the underestimation of error covariance is avoided, and the estimation accuracy is improved.
Overall, the results in Figure 6 provide the ADCS performance guideline during the system design review stage. The expected quaternion estimation error increment from Figure 6 allows the evaluation if an additional sensor, such as a star tracker is required to meet the system design requirement. STAR centre is presently developing a 50 kg microsatellite with a star tracker for high-precision pointing applications. The applicability of SBPKF will be verified on the engineering model. Subsequently, its in-orbit attitude determination accuracy performance will be benchmarked with the quaternion computed by the star tracker. Table 3 compares the number of multiplications required by EKF, SBPKF, MAEKF, and MC-KF at each iteration. The number of multiplication required for matrix multiplication and an inverse matrix is based on the method in [20]. The MAEKF's multiplication number is the average multiplication number per sampling step within one second. This is due to the reason that the number of available measurement pairs for MAEKF varies between 1 and 2 at each sampling step. For the SBPKF, MAEKF, and MC-KF, we assume the scenario where δt > 5 ms. Table 3 shows that the number of multiplications required by MC-KF and MAEKF is approximately two times higher than EKF. The number of multiplication required by SBPKF is 17% lesser than MAEKF and 29% lesser than MC-KF, respectively. Although the SBPKF requires 65% more multiplication compared to EKF, the quaternion estimation error is improved by 30 to 45% with a similar gyro bias estimation error. Thus, the SBPKF achieves a higher attitude estimation accuracy at an expense of higher computational complexity.

Conclusions
This paper presented the simplex-back-propagation Kalman filter (SBPKF) for delayed sensor data-based quaternion and gyro bias estimation. The estimated measurements are formulated with the consideration of sensor data's time delay to improve the quaternion estimation accuracy. In addition, the noise covariance matrices are derived with the presence of navigation error, gyro bias, and gyro noises to prevent underestimating the estimation error. The detailed implementation of the moving-covariance Kalman filter (MC-KF) for quaternion estimation has been presented. Monte Carlo simulations have been conducted to compare the accuracy and computational complexity of the proposed SPBKF, extended Kalman filter (EKF), modified adaptive EKF (MAEKF), and MC-KF. The results show that the SBPKF average estimated quaternion error is at least 30% lower than both EKF and MC-KF. Although the MAEKF quaternion estimation error is slightly higher than the SBPKF, its gyroscope bias estimation error is three times higher than SBPKF. In addition, the SPBKF computational complexity is 17% better than MAEKF and 29% better than MC-KF.
For future work, the SPBKF will be evaluated using the engineering model of a 50 kg microsatellite. The impact of the digital signal processor's 8-bit floating point on the overall SPBKF accuracy performance will also be investigated.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
This appendix provides an overview of the derivation of the Earth's magnetic field and its associate partial derivative. From Appendix H in [38], the North-East-Down components of the Earth's magnetic field vector (or B θ , B φ and B r ) are given as Here, θ c is the geocentric co-latitude such that where θ is the geodetic latitude. The partial derivate of each component in B, with respect to the North-East-Down direction (or known as θ, φ, and r, respectively) are given as It is noted that based on (A2), the partial derivate of the Earth's magnetic field with respect to geodetic latitude θ in (A3) to (A5) has been simplified based on the following definition ∂B ∂θ In addition, from Appendix H in [38], the Gauss function P n,m (θ c ) are given as below: P 0,0 (θ c ) = 1 (A7) P n,n (θ c ) = sin θ c P n−1,n−1 (θ c ) (A8) P n,m (θ c ) = cos θ c P n−1,m (θ c ) − K n,m P P n−2,m (θ c ) (A9) K n,m P ≡ (n−1) 2 −m 2 (2n−1)(2n−3) if n > 1 0 if n = 1 (A10) The first and second order partial derivative with respect to θ c can be obtained via (A7) to (A9).