Satellite Inertia Parameters Estimation Based on Extended Kalman Filter

ABSTRACT:

The moment of inertia parameters play a critical role in assuring the spacecraft mission throughout its lifetime. However, determination of the moment of inertia is a key challenge in operating satellites. During satellite mission, those parameters can change in orbit for many reasons such as sloshing, fuel consumption, etc. Therefore, the inertia matrix should be estimated in orbit to enhance the attitude estimation and control accuracy. This paper investigates the use of gyroscope to estimate the attitude rate and inertia matrix for low earth orbit satellite via extended Kalman filter. Simulation results show the effectiveness and advantages of the proposed algorithm in estimating these parameters without knowing the nominal inertia. The robustness of the proposed algorithm has been validated using the Monte-Carlo method. The obtained results demonstrate that the accuracy of the estimated inertia and angular velocity parameters is satisfactory for satellite with coarse accuracy mission requirements. The proposed method can be used for different types of satellites.

KEYWORDS:
Satellite inertia tensor estimation; Extended Kalman filter; Spacecraft inertia parameters; Gyroscope

INTRODUCTION

Spacecraft inertia parameters are prone to unpredictable changes throughout a mission due to many reasons such as fuel consumption, fuel sloshing, deployment failures, collision with unexpected object, damage, connection with other parts, deorbiting, structure changes or any events which substantially influence the inertia matrix (Kim et al. 2010; 2016; Manchester and Peck 2017). Knowledge of the correct inertia parameters of a spacecraft is vital to its control system.

In-orbit estimation of the spacecraft’s inertia parameters is more efficient than calculating these parameters on-ground before launch because the precision of the latter is limited. Several researches and engineering projects have been proposed to estimate full components of inertia matrix, which consist of moment of inertia (MOI) and product of inertia (POI) elements.

Least squares algorithms have been considered in various research works for identifying satellite inertia parameters. Tanygin and Williams (1997) utilize a standard least square; this method allows determination of both the inertia matrix and the position of the center of mass. In Bordany et al. (2000), a recursive least square (RLS) procedure is proposed to identify satellite inertia matrix and thruster parameters in orbit for mini-satellite UoSAT-12 built by Surrey Satellite Technology Ltd. In Psiaki (2005), an algorithm is developed using the inner and outer least-squares optimizations to estimate the dynamic model parameters. Keim et al. (2006) have presented an inertia estimation technique based on a constraint least squares minimization problem with LMIs (linear matrix inequalities). Most of these algorithms require integrals or derivatives of measured spacecraft state variables as input, necessitating prefiltering that introduces additional complexity and numerical error.

Many researches using Kalman filter for estimating spacecraft’s inertia parameters were proposed in the literature. However, extended Kalman filter (EKF) is more suitable for nonlinear systems. Kutlu et al. (2007) have combined Kalman Filter algorithms for attitude determination part using gyroscopes, magnetometer and a sun sensor, and for parameter identification part using a Least Square estimation algorithm. Schwartz and Hall (2003) have presented a comparison and analyses for estimating simultaneously the states and parameters using EKF and unscented Kalman Filter (UKF). In (Kim et al. 2010) the filtered data of angular rates are obtained using the EKF and they are used to the batch method for estimating inertia properties of STSAT-3; then Savitzky-Golay filter (SGF) was proposed in (Kim et al. 2014) to estimate full inertia properties and recently this filter was combined with the extended Kalman filter (EKF) using a gyroscope measurements (Kim et al. 2016). In Inamori et al. (2011) an EKF is proposed for the estimation and compensation of a low frequency residual magnetic moment of micro-astronomy satellite Nano-JASMINE developed at the University of Tokyo. This method was recently used by Söken and Hajiyev (2014) for estimation of pico-satellite attitude dynamics and external torques via unscented Kalman filter.

The aim of this paper is to estimate the spacecraft’s inertia components, MOI and POI, of the inertia matrix based on EKF using the algorithm proposed in (Inamori et al. 2011; Söken and Hajiyev 2014). The proposed estimator does not use any nominal inertia matrix and the initialization will be random. The advantage of not using any nominal inertia matrix is to reduce the impact of the error of the inertia matrix on the angular velocity estimation. In addition, the noise, the bias caused by gyroscope and the external disturbance torque vector are taken in consideration. Monte Carlo method is used to test the effectiveness and robustness of the proposed estimator.

The remainder of the paper is organized as follows: the proposed method and the different steps of the EKF design are described in section Estimation Method; section Simulation and Results presents the simulation results with discussions; finally, the last section concludes this paper.

ESTIMATION METHOD

State Vector and System Model

The state vector to be estimated, x, which has nine elements is expressed as (Eq. 1)

(1) $\mathit{x}={\left[{\mathit{\omega }}_{\mathit{x}}{\mathit{\omega }}_{\mathit{y}}{\mathit{\omega }}_{\mathit{z}}{\mathit{I}}_{\mathit{xx}}{\mathit{I}}_{\mathit{yy}}{\mathit{I}}_{\mathit{zz}}{\mathit{I}}_{\mathit{xy}}{\mathit{I}}_{\mathit{xz}}{\mathit{I}}_{\mathit{yz}}\right]}^{\mathit{T}}$

where ω = [ω_x ω_y ω_z]^T; J_m = [I_xx I_yy I_zz]^T; J_p = [I_xy I_xz I_yz]^T. ω is the angular velocity vector in the inertial frame; J_m and J_p are the moment of inertia (MOI) and product of inertia (POI) vector respectively.

The dynamics of the spacecraft in inertial space will be governed by Euler’s equations of motion (Wertz 1978). With the added influence of the external disturbance torque, applied control torque and since we estimate the constant components of inertia matrix, the dynamics model of the spacecraft can be expressed as follows in vector form as (Eq. 2),

(2) $\left[\begin{array}{c}\dot{\mathit{\omega }}\hfill \\ {\dot{\mathit{J}}}_{\mathit{m}}\hfill \\ {\dot{\mathit{J}}}_{\mathit{p}}\hfill \end{array}\right]=\left[\begin{array}{c}{\mathit{J}}^{-1}\left[{\mathit{T}}_{\mathit{c}}+{\mathit{T}}_{\mathit{d}}-\left(\mathit{\omega }\times \mathit{J}\mathit{\omega }\right)\right.\\ -\frac{1}{{\mathit{\tau }}_{\mathit{m}}}{\left[\mathit{I}\right]}_{3\times 3}{\mathit{J}}_{\mathit{m}}\\ -\frac{1}{{\mathit{\tau }}_{\mathit{p}}}{\left[\mathit{I}\right]}_{3\times 3}{\mathit{J}}_{\mathit{p}}\end{array}\right]$

where J is the inertia matrix of satellite; T_c and T_d are the torque command and external disturbance torque vector respectively; I is the identity matrix; τ_m and τ_p are the time constant of the MOI and POI vector, respectively.

The full non-linear model and the state equation can be obtained using Eq. 3,

(3) $\dot{\mathit{x}}\left(\mathit{t}\right)=\mathit{f}\left\{\mathit{x}\left(\mathit{t}\right),\mathit{t}\right\}+\mathit{w}\left(\mathit{t}\right)$

where w(t) is zero mean white system noise with covariance matrix Q (Eq. 4)

(4) $\mathit{w}\left(\mathit{t}\right)=\mathit{N}\left\{0,\mathit{Q}\left(\mathit{t}\right)\right\}$

State Transition Matrix

A discrete version of the EKF will be used to enable onboard implementation, for short sampling period Δt the discrete state transition matrix Φ can be approximated using a first order Taylor series expansion as (Eq. 5) (Gelb 1989; Brown and Hwang 1997),

(5) ${\mathit{\Phi }}_{\mathit{k}}\approx {\mathit{I}}_{9\times 9}+\left[\mathit{F}\left\{\hat{\mathit{x}}\left({\mathit{t}}_{\mathit{k}}\right),{\mathit{t}}_{\mathit{k}}\right\}\Delta \mathit{t}\right]$

The state matrix $\mathit{F}\left\{\hat{\mathit{x}}\left({\mathit{t}}_{\mathit{k}}\right),{\mathit{t}}_{\mathit{k}}\right\}$ is obtained by partially differentiating the dynamics (Eq. 2) with respect to the state vector (Eq. 1). $\mathit{F}\left\{\hat{\mathit{x}}\left({\mathit{t}}_{\mathit{k}}\right),{\mathit{t}}_{\mathit{k}}\right\}$ is expressed in the matrix form as (Eqs. 6 and 7),

(6) $\mathit{F}\left\{\hat{\mathit{x}}\left({\mathit{t}}_{\mathit{k}}\right),{\mathit{t}}_{\mathit{k}}\right\}=\left[\begin{array}{ccc}\frac{\partial \dot{\mathit{\omega }}}{\partial \mathit{\omega }}& \frac{\partial \dot{\mathit{\omega }}}{\partial {\mathit{J}}_{\mathit{m}}}& \frac{\partial \dot{\mathit{\omega }}}{\partial {\mathit{J}}_{\mathit{p}}}\\ \frac{\partial {\dot{\mathit{J}}}_{\mathit{m}}}{\partial \mathit{\omega }}& \frac{\partial \dot{\mathit{J}}}{\partial {\mathit{J}}_{\mathit{m}}}& \frac{\partial {\dot{\mathit{J}}}_{\mathit{m}}}{\partial {\mathit{J}}_{\mathit{p}}}\\ \frac{\partial {\dot{\mathit{J}}}_{\mathit{p}}}{\partial \mathit{\omega }}& \frac{\partial {\dot{\mathit{J}}}_{\mathit{p}}}{\partial {\mathit{J}}_{\mathit{m}}}& \frac{\partial {\dot{\mathit{J}}}_{\mathit{p}}}{\partial {\mathit{J}}_{\mathit{p}}}\end{array}\right]$

(7) $\mathit{F}\left\{\hat{\mathit{x}}\left({\mathit{t}}_{\mathit{k}}\right),{\mathit{t}}_{\mathit{k}}\right\}=\left[\begin{array}{ccc}\frac{\partial \dot{\mathit{\omega }}}{\partial \mathit{\omega }}& \frac{\partial \dot{\mathit{\omega }}}{\partial {\mathit{J}}_{\mathit{m}}}& \frac{\partial \dot{\mathit{\omega }}}{\partial {\mathit{J}}_{\mathit{p}}}\\ 0_{3\times 3}& -\frac{1}{{\mathit{\tau }}_{\mathit{m}}}{\left[\mathit{I}\right]}_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& -\frac{1}{{\mathit{\tau }}_{\mathit{p}}}{\left[\mathit{I}\right]}_{3\times 3}\end{array}\right]$

For convenience, let the matrix $\mathit{F}\left\{\hat{\mathit{x}}\left({\mathit{t}}_{\mathit{k}}\right),{\mathit{t}}_{\mathit{k}}\right\}$ be (Eq. 8)

(8) $\mathit{F}\left\{\hat{\mathit{x}}\left({\mathit{t}}_{\mathit{k}}\right),{\mathit{t}}_{\mathit{k}}\right\}=\left[\begin{array}{ccc}{\mathit{F}}_{3\times 3}^{\mathit{\omega }\mathit{\omega }}& {\mathit{F}}_{3\times 3}^{\mathit{\omega }{\mathit{J}}_{\mathit{m}}}& {\mathit{F}}_{3\times 3}^{\mathit{\omega }{\mathit{J}}_{\mathit{r}}}\\ 0_{3\times 3}& {\mathit{F}}_{3\times 3}^{{\mathit{J}}_{\mathit{m}}{\mathit{J}}_{\mathit{m}}}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& {\mathit{F}}_{3\times 3}^{{\mathit{J}}_{\mathit{r}}{\mathit{J}}_{\mathit{r}}}\end{array}\right]$

Process Noise Covariance Matrix

The process noise covariance matrix Q is defined as (Eq. 9) (Gelb 1989; Brown and Hwang 1997)

(9) ${\mathit{Q}}_{\mathit{k}}=\underset{0}{\overset{\Delta \mathit{t}}{\int }}\underset{0}{\overset{\Delta \mathit{t}}{\int }}\mathit{\Phi }\left(\Delta \mathit{t},\mathit{u}\right)\mathit{E}\left[\mathit{w}\left(\mathit{u}\right)\mathit{w}{\left(\mathit{v}\right)}^{\mathit{T}}\right]\mathit{\Phi }{\left(\Delta \mathit{t},\mathit{v}\right)}^{\mathit{T}}\mathit{dudv}$

where u and v are independent variables; for our application, it is assumed that only the angular velocity terms have process noise, and the matrix E of dimension (9 × 9) is given by (Eq. 10):

(10) $\mathit{E}\left[\mathit{w}\left(\mathit{u}\right)\mathit{w}{\left(\mathit{v}\right)}^{\mathit{T}}\right]=\left[\begin{array}{cc}{\mathit{E}}_{11}& 0_{3\times 6}\\ 0_{6\times 3}& 0_{6\times 6}\end{array}\right]$

where (Eq. 11),

(11) ${\mathit{E}}_{11}=\mathit{diag}\left[{\mathit{s}}_{\mathit{x}}\mathit{\delta }\left(\mathit{u}-\mathit{v}\right){\mathit{s}}_{\mathit{y}}\mathit{\delta }\left(\mathit{u}-\mathit{v}\right){\mathit{s}}_{\mathit{z}}\mathit{\delta }\left(\mathit{u}-\mathit{v}\right)\right]$

δ is the Dirac Delta function and s is the spectral amplitude for the angular velocity.

In Eq. 9 the term Φ (Δt, u)E[w(u) w(v)^T] Φ (Δt, v)^T is computed as (Eq. 12)

(12) $\left[\begin{array}{ccc}{\mathit{I}}_{3\times 3}+{\mathit{F}}_{3\times 3}^{\mathit{\omega }\mathit{\omega }}\mathit{u}& {\mathit{F}}_{3\times 3}^{\mathit{\omega }{\mathit{J}}_{\mathit{m}}}\mathit{u}& {\mathit{F}}_{3\times 3}^{\mathit{\omega }{\mathit{J}}_{\mathit{r}}}\mathit{u}\\ 0_{3\times 3}& {\mathit{I}}_{3\times 3}+{\mathit{F}}_{3\times 3}^{{\mathit{J}}_{\mathit{m}}{\mathit{J}}_{\mathit{m}}}\mathit{u}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& {\mathit{I}}_{3\times 3}+{\mathit{F}}_{3\times 3}^{{\mathit{J}}_{\mathit{r}}{\mathit{J}}_{\mathit{r}}}\mathit{u}\end{array}\right]\times \left[\begin{array}{ccc}{\mathit{E}}_{11}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\end{array}\right]\times \left[\begin{array}{ccc}{\mathit{I}}_{3\times 3}+{\mathit{F}}_{3\times 3}^{\mathit{\omega }\mathit{\omega }}\mathit{v}& 0_{3\times 3}& 0_{3\times 3}\\ {\mathit{F}}_{3\times 3}^{\mathit{\omega }{\mathit{J}}_{\mathit{m}}}\mathit{v}& {\mathit{I}}_{3\times 3}+{\mathit{F}}_{3\times 3}^{{\mathit{J}}_{\mathit{m}}{\mathit{J}}_{\mathit{m}}}\mathit{v}& 0_{3\times 3}\\ {\mathit{F}}_{3\times 3}^{\mathit{\omega }{\mathit{J}}_{\mathit{r}}}\mathit{v}& 0_{3\times 3}& {\mathit{I}}_{3\times 3}+{\mathit{F}}_{3\times 3}^{{\mathit{J}}_{\mathit{r}}{\mathit{J}}_{\mathit{r}}}\mathit{v}\end{array}\right]$

Using definition and proprieties of Dirac Delta function; the resulted process noise covariance matrix is expressed as (Eq. 13)

(13) ${\mathit{Q}}_{\mathit{k}}=\left[\begin{array}{cc}{\mathit{q}}_{3\times 3}& 0_{3\times 6}\\ 0_{6\times 3}& 0_{6\times 6}\end{array}\right]$

where (Eq. 14)

(14) $\begin{array}{c}{\mathit{q}}_{3\times 3}={\mathit{s}}_{\mathit{\omega }}\Delta \mathit{t}+{\mathit{s}}_{\mathit{\omega }}{\left({\mathit{F}}_{3\times 3}^{\mathit{\omega }\mathit{\omega }}\right)}^{\mathit{T}}\frac{{\left(\Delta \mathit{t}\right)}^2}{2}+{\mathit{s}}_{\mathit{\omega }}{\mathit{F}}_{3\times 3}^{\mathit{\omega }\mathit{\omega }}\frac{{\left(\Delta \mathit{t}\right)}^2}{2}+{\mathit{F}}_{3\times 3}^{\mathit{\omega }\mathit{\omega }}{\mathit{s}}_{\mathit{\omega }}{\left({\mathit{F}}_{3\times 3}^{\mathit{\omega }\mathit{\omega }}\right)}^{\mathit{T}}\frac{{\left(\Delta \mathit{t}\right)}^3}{3}\\ {\mathit{s}}_{\mathit{\omega }}=\left[\begin{array}{ccc}{\mathit{s}}_{\mathit{x}}& 0& 0\\ 0& {\mathit{s}}_{\mathit{y}}& 0\\ 0& 0& {\mathit{s}}_{\mathit{z}}\end{array}\right]\end{array}$

Gyroscope Observation Equation

The EKF use the angular rate measurement from the gyroscope and the observation model is given as follows (Eq. 15):

(15) ${\mathit{z}}_{\mathit{k}}={\mathit{\omega }}_{\mathit{k}}+{\mathit{b}}_{\mathit{k}}+{\mathit{\eta }}_{\mathit{k}}$

where (Eq. 16)

(16) ${\mathit{\eta }}_{\mathit{k}}=\mathit{N}\left\{0,{\mathit{R}}_{\mathit{k}}\right\}$

η_k is the Gaussian measurement noise with zero mean and a covariance matrix R; b denotes biases in each axis that is a slowly variable function over time. Therefore the gyroscope observation matrix H with dimension (3 × 9) becomes (Eq. 17):

(17) ${\mathit{H}}_{\mathit{k}}=\left[\frac{\partial \mathit{z}}{\partial \mathit{\omega }}\frac{\partial \mathit{z}}{\partial {\mathit{J}}_{\mathit{m}}}\frac{\partial \mathit{z}}{\partial {\mathit{J}}_{\mathit{r}}}\right]=\left[{\mathit{I}}_{3\times 3}{0}_{3\times 6}\right]$

EKF Algorithm

The estimation of the EKF operates in two primary cycles, propagation and correction. During the propagation cycle, the filter propagates the state of the system, using a system model to predict the state of the system one-time step in the future. The correction cycle inputs measurements of the system state and utilizes these observations to correct for differences between the state propagated from the system model and the measured satellite state.

Propagation Cycle:

• 1. Propagation of the state equation (Eq. 2), by using the second Adams integration method (Eq. 18)

  (18) ${\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}={\hat{\mathit{x}}}_{\mathit{k}-1/\mathit{k}-1}+\frac{1}{2}\left(3{\dot{\mathit{x}}}_{\mathit{k}-1/\mathit{k}-1}-{\dot{\mathit{x}}}_{\mathit{k}-2/\mathit{k}-2}\right)\Delta \mathit{t}$

where Δt is the integration step size.

• 2. Propagation of the covariance estimate matrix using the state transition matrix (Eq. 19)

  (19) ${\mathit{P}}_{\mathit{k}/\mathit{k}-1}={\mathit{\Phi }}_{\mathit{k}-1}{\mathit{P}}_{\mathit{k}-1/\mathit{k}-1}{\mathit{\Phi }}_{\mathit{k}-1}^{\mathit{T}}+{\mathit{Q}}_{\mathit{k}-1}$

Correction Cycle:

1. Compute Observation Matrix (H) and update Kalman Gain (K) (Eq. 20)

   (20) ${\mathit{K}}_{\mathit{k}}={\mathit{P}}_{\mathit{k}/\mathit{k}-1}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}{\left[{\mathit{H}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}/\mathit{k}-1}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}+{\mathit{R}}_{\mathit{k}}\right]}^{-1}$

2. Update state estimate with Measurement z_k (Eq. 21)

   (21) ${\mathit{x}}_{\mathit{k}/\mathit{k}}={\mathit{x}}_{\mathit{k}/\mathit{k}-1}+{\mathit{K}}_{\mathit{k}}\left[{\mathit{z}}_{\mathit{k}}-{\mathit{H}}_{\mathit{k}}{\mathit{x}}_{\mathit{k}/\mathit{k}-1}\right]$

3. Update covariance estimate matrix (Eq. 22)

   (22) ${\mathit{P}}_{\mathit{k}/\mathit{k}}=\left[{\mathit{I}}_{9\times 9}-{\mathit{K}}_{\mathit{k}}{\mathit{H}}_{\mathit{k}}\right]{\mathit{P}}_{\mathit{k}/\mathit{k}-1}$

Observability

The observability is a structural characteristic of a state representation of a system that indicates the ability of a system to determine the history of a state from knowledge of the output variables measured. The notion of observability involves the dynamic and the output matrix; the observability criterion (Kalman) is based on the rank test of the observability matrix OB, which is defined as (Eq. 23) (Simon 2006),

(23) $\mathit{OB}={\left[\mathit{H}\mathit{HF}{\mathit{HF}}^2{\mathit{HF}}^3{\mathit{HF}}^4{\mathit{HF}}^5{\mathit{HF}}^6{\mathit{HF}}^7{\mathit{HF}}^8\right]}^{\mathit{T}}$

The estimation system is observable if the rank of the matrix OB is equal to size state vector. Figure 1 shows that the rank of the observability matrix is equal to 9. Therefore, the estimation system is observable.

Figure 1
Rank of observability matrix.

Extended Kalman and Savitzky-Golay Filters

In Kim et al. (2016), a combined methodology is proposed to estimate full moment of inertia including the off diagonal of inertia matrix. The estimation process consists of the following three steps: first, the noise measurement is filtered using the EKF; next, the reliable angular acceleration is calculated using the Savitzky-Golay filter (SGF) for an even number of sampled data; last, the moment of inertia matrix is estimated using the linear least squares (LLS) based on the suggested regressor matrix. In this paper, for simplification this estimator is named EKF#1.

The main differences between the proposed EKF and EKF#1 will be highlighted. The covariance and the true state are predicted in the propagation cycle. However, all matrices involved in the covariance propagation are (9 × 9) dimension, using EKF, and are (3 × 3) dimension using EKF#1, which significantly reduces the computational load upon this filter. In the proposed EKF, the state vector is propagated and directly updated. However, the EKF#1 process is more complex compared to the proposed EKF because it uses the SGF to calculate the angular accelerations, linear least squares (LLS) and the regressor matrix to estimate the inertia vector which increases the computational load. In addition, the EKF#1 uses nominal inertia matrix in propagation cycle, the estimated value of the angular velocity is always calculated using this nominal matrix. In the proposed EKF, the angular velocity is calculated using the estimated inertia, which could reduce the estimated angular velocity error.

SIMULATION AND RESULTS

Our simulation is based on an attitude rate propagator for low earth orbit microsatellite, which models the actual satellite dynamics and outputs Euler rates over a specified period of time. The simulation program is developed from dynamic model to propagate expected rates. Due to the construction of this model, it is quite easy to simplify the model to take account of all disturbances torques. An extended Kalman filter (EKF) is used to estimate the nine-state vector from gyroscope measurement. The nine-state EKF vector is composed of three-element body rates vector combined with six-element inertia vector. The simulation was run for 600 s and a summary of the operating constants is presented in Table 1, whereas the EKF parameters are presented in Table 2.

Table 1
Simulation parameters.

Parameters		Value
Orbital parameters	Inclination (deg)	97.8
	Initial ω (deg/s)	[0 -0.06 0]T
	Altitude (km)	600
Sample time (s)		10
True MOI (kg∙m2)		[14.2 17.3 20.3]T
True POI (kg∙m2)		[0.0867 0.1357 0.6016]T
Command torque (Nm)		[1 1 -2]T × 10-3

Table 2
EKF parameters.

Parameters	Value
Initial ω (deg/s)	[0 -2.3 1.15]T
Initial MOI (kg∙m2)	[25 25 25]T
Initial POI (kg∙m2)	[2 2 2]T
Gyroscope noise level (deg/s)	0.0104
[τm τp]	[106 105]

Figure 2 presents the estimated angular velocity. The filtered angular velocities and 3σ boundaries are shown in Fig. 3. The angular velocity error result indicates that the global root mean square error (RMSE) is about 3.6 × 10^-3 deg/s. Figure 4 illustrates the errors variances of the estimated angular velocity in three-axis. The results show that the standard deviation is approximately 2.4 × 10^-3 deg/s.

Figure 2
Estimated angular velocity.

Figure 3
Estimated angular velocity error

Figure 4
Error variance of the estimate angular velocity.

Figures 5 and 6 show the estimated MOI tensor and POI tensor respectively; it can also be seen that the convergence time of estimated inertia vector on the true inertia vector is about 350 s. Table 3 summarizes the results of this estimation and comparison with estimator EKF#1 (Kim et al. 2016). Both filters proved to meet the requirements ADCS inertia tensor estimation success with little difference in the accuracy obtained. EKF#1 has a better accuracy than the proposed EKF for the POI estimate, but for the estimation of MOI the proposed EKF is better. In addition, we didn’t use any nominal inertia matrix and the difference between the initiate and the true inertia matrix is more than 50% and the proposed EKF converges quickly as EKF#1.

Figure 5
Estimated MOI tensor.

Figure 6
Estimated POI tensor.

Table 3
EKF estimation results.

Inertia	True (kg∙m2)	EKF		EKF#1
		Estimated (kg∙m2)	Error (%)	Estimated (kg∙m2)	Error (%)
Ixx	14.20	14.1905	0.06	14.1823	0.12
Iyy	17.30	17.2954	0.02	17.2876	0.07
Izz	20.30	20.3017	0.008	20.2859	0.07
Ixy	0.0867	0.0871	0.46	0.0864	0.39
Ixz	0.1357	0.1345	0.88	0.1350	0.48
Iyz	0.6016	0.6001	0.25	0.6009	0.12

The accuracy performance and robustness of the proposed algorithm are validated using Monte-Carlo method, which consists of 10000 Monte-Carlo runs. This simulation is divided into five steps; each step contains 2000 runs for an attenuation of the inertia matrix of 10%, the last step the inertia matrix is decreased to 50%. For each Monte-Carlo run, the initial angular velocity and initial inertia tensor matrix of estimator was picked randomly; the starting angular velocity, MOI and POI came from a uniform population between ±0.06 deg·s^-1, ±30 kg·m² and ±2 kg·m², respectively. The accuracy characteristics of the algorithm are captured in Figs. 7 and 8, which display a histogram of MOI and POI error (%) for the 10000 Monte Carlo running respectively. The algorithm never diverged in all Monte-Carlo runs. Figures 7 and 8 confirm the estimation capability of the proposed estimator with regard to the variation of the most critical parameters.

Figure 7
Histogram of MOI error for 10000 Monte-Carlo runs.

Figure 8
Histogram of POI error for 10000 Monte-Carlo runs.

CONCLUSION

In this paper, we have presented an algorithm to estimate the moment of inertia, product of inertia elements and angular velocity based on extended Kalman filter using only gyroscope measurement applied to a LEO microsatellite. The proposed estimation algorithm is designed and compared to the batch least squares method combined with the EKF. The simulation results show that the proposed algorithm based on EKF is faster and has a better accuracy for the estimation of moment of inertia without using any nominal inertia matrix. In addition, Monte-Carlo method has been used with different scenarios to validate the proposed algorithm. The simulation results demonstrate the effectiveness and robustness of the presented algorithm without using any information of the nominal inertia. Indeed, the proposed filter can guarantee the estimation performance up to 50% error of the true values of the inertia matrix. For future work, we aim to develop a new filter algorithm that estimates simultaneously the full inertia tensor and disturbance torque.

• Bellar A https://orcid.org/0000-0003-2786-7627
• Si Mohammed MA https://orcid.org/0000-0001-6561-4757

• FUNDING
  There are no funders to report.

REFERENCES

• Bordany RE, Steyn WH, Crawford M (2000) In-orbit estimation of the inertia matrix and thruster parameters of UoSAT-12. Presented at: 14th Annual AIAA/USU Conference on Small Satellites; Logan, USA.
• Brown RG, Hwang PYC (1997) Introduction to random signals and applied Kalman filtering. New York: John Willey.
• Gelb A (1989) Applied optimal estimation. Cambridge: MIT Press.
• Inamori T, Sako N, Nakasuka S (2011) Compensation of time-variable magnetic moments for a precise attitude control in nano- and micro-satellite missions. Advances in Space Research 48(3):432-440. https://doi.org/10.1016/j.asr.2011.03.036
  » https://doi.org/10.1016/j.asr.2011.03.036
• Keim JA, Acikmese AB, Shields JF (2006) Spacecraft inertia estimation via constrained least squares. Presented at: IEEE Aerospace Conference; Big Sky, USA. https://doi.org/10.1109/AERO.2006.1655995
  » https://doi.org/10.1109/AERO.2006.1655995
• Kim D, Lee S, Leeghim H (2014) Inertia estimation using Savitzky-Golay filter. Presented at: KSAS Fall Conference; Jeju, Republic of Korea.
• Kim D, Yang S, Lee S (2016) Rigid body inertia estimation using extended Kalman and Savitzky-Golay filters. Mathematical Problems in Engineering 2016:2962671. https://doi.org/10.1155/2016/2962671
  » https://doi.org/10.1155/2016/2962671
• Kim D, Yang S, Cheon DI, Lee S, Oh HS (2010) Combined estimation method for inertia properties of STSAT-3. Journal of Mechanical Science and Technology 24(8):1737-1741. https://doi.org/10.1007/s12206-010-0521-2
  » https://doi.org/10.1007/s12206-010-0521-2
• Kutlu, A., Haciyev, C. and Tekinalp, O., 2007 Attitude determination and rotational motion parameters identification of a LEO satellite through magnetometer and sun sensor data. Presented at: 3rd International Conference on Recent Advances in Space Technologies; Istanbul, Turkey. https://doi.org/10.1109/RAST.2007.4284033
  » https://doi.org/10.1109/RAST.2007.4284033
• Manchester ZR, Peck MA (2017) Recursive inertia estimation with semidefinite programming. Presented at: AIAA Guidance, Navigation, and Control Conference; Gaylord, USA. https://doi.org/10.2514/6.2017-1902
  » https://doi.org/10.2514/6.2017-1902
• Psiaki ML (2005) Estimation of a spacecraft’s attitude dynamics parameters by using flight data. Journal of Guidance, Control, and Dynamics 28(4):594-603. https://doi.org/10.2514/1.7362
  » https://doi.org/10.2514/1.7362
• Schwartz JL, Hall CD (2003) Comparison of system identification techniques for a spherical air-bearing spacecraft simulator. Presented at: AAS/AIAA Astrodynamics Specialist Conference; Big Sky, USA.
• Simon D (2006) Optimal state estimation: Kalman, H∞, and nonlinear approaches. Chichester: John Wiley & Sons.
• Söken HE, Hajiyev C (2014) Estimation of pico-satellite attitude dynamics and external torques via unscented Kalman filter. Journal of Aerospace Technology and Management 6(2):149-157. https://doi.org/10.5028/jatm.v6i2.352
  » https://doi.org/10.5028/jatm.v6i2.352
• Tanygin S, Williams T (1997) Mass Property Estimation Using Coasting Maneuvers. Journal of Guidance, Control, and Dynamics 20(4):625-632. https://doi.org/10.2514/2.4099
  » https://doi.org/10.2514/2.4099
• Wertz JR (1978) Spacecraft Attitude Determination and Control. Dordrecht: Springer.

Edited by

Section Editor: Luiz Martins-Filho

Publication Dates

• Publication in this collection
  28 Mar 2019
• Date of issue
  2019

History

• Received
  12 June 2017
• Accepted
  12 June 2018