# Modeling and Prediction of Momentum Wheel Speed Data

Authors: Jichao Li; Xiaxia Wang; Chaobo Chen; Song Gao
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

To solve the problems of data loss and unequal interval of momentum wheel (MW) speed during a satellite stable operation, this paper presents a multidimensional AR model. A Lagrange interpolation method is used to convert measurements to equal interval data, and the FFT algorithm is adopted to calculate the period of MW speed variation. The long data sequence is converted into multidimensional time series, based on the equal interval data and the period. A multidimensional AR model is established, and the least square method is used to estimate the model parameters. The future data trend is predicted by the proposed model. Simulation results show that the prediction algorithm can achieve the across cycle prediction of the MW speed data.


---

## Body

## 1. Introduction

The momentum wheel (MW) is the key actuator of the Satellite Attitude Control System (SACS), which includes bearing assembly, motor, casing, and wheels. In the long running of the three-axis stabilized GEO satellite, the MW speed data has a periodical change characteristic. Through long-term monitoring for MW speed of a three-axis stabilized GEO satellite, [1] found that its variation can be sorted into the short-period terms (in days) and the long-period terms (in years), and the periodic variation of the MV is also analyzed. The equation of the MW speed is obtained by analyzing the operation principle of the bias MW, and the graphic of MW speed change is described in [2]. In [3], the authors considered that the reliability information of MW can be provided by collecting degradation data when there exists certain performance characteristics that degrade over time; they develop a reliability modeling and life estimation approach for MW used in satellites based on the expectation maximization (EM) algorithm from a Wiener degradation model. GmbH has installed 588 MWs in 235 satellites, which is required fault-free operation more than 2200 years [4].Under the complex working environment, such as water mist, salt spray, corrosion, and aging resulting from fuels containing sodium, sulfur, vanadium, etc. which are burned, the dynamic characteristics of speed sensors are severely influenced, which causes the problems that the data detected is inaccurate or could not be detected in a certain time that means false detection or data loss [5]. Furthermore, the data detected may be equal interval incompletely because of the unstable precision of a sensor. As we know, the rotational speed is an important state parameter of the momentum wheel. Once the speed sensors occur to faults, it will cause a wrong operation of control system and even lead to accidents. The Telstar satellite of America was lost early on 21 February 1963 because of the abnormal data of a sensor caused by cosmic rays. This was the first spacecraft loss due to radiation effects [6].The remainder of this paper is organized as follows. In Section2, we briefly analyze the change characteristic of the MW speed data, using Lagrange interpolation and Fast Fourier Transform Algorithm (FFT) to determine the period of MW speed data. Section 3 establishes the multidimensional model to predict the MW speed data and verifies the validity of the model. Finally, the paper concludes with some brief, summary remarks for the present work.


## 2. Obtainment of the Variation Period of the Momentum Wheel Date

### 2.1. Analysis of the Momentum Wheel Speed Data

The momentum wheel speed data from a simulation test of communication satellite is time series data. Assuming that the observed MW speed isyti at time ti, written as yi, the change curve of MW speed is shown in Figure 1. Because of accidental factors such as cosmic ray and accidental contact failure of circuit, the sensor data is abnormal [7]; for example, the observed data is incompletely equal interval distribution and the partially observed data loss. Meanwhile, it can be seen that the MW speed is approximately periodic and its amplitude grows linearly. In order to predict MW speed, the Lagrange interpolation method [8] is used to calculate the MW speed, so that the unequal interval data is converted into equal interval data.Figure 1
The change curves of momentum wheel speed.


### 2.2. Preprocessing of Equal Interval Interpolation of Momentum Wheel Speed Data

The time series data of the observed MW speed isY=t1,y1,t2,y2,⋯,tN,yNT, where N is the total number of data. However, Y is the unequal interval time series data, and the forecast model based on equal interval time series data. So Y has to convert into equal interval time series data by the interpolation method.Firstly, set the interpolation intervalp. To interpolate the data segment from t1,y1 to tm+1,ym+1 in Y, the Lagrange interpolation formula [9] can be given by
(1)Lmt=∑k=1m+1∏i=0i≠0m+1t−titk−tiyk,where t is the interpolation time; t1,t2,⋯,tm+1 are the observation time which satisfy t1<t2<⋯<tm+1.According to formula (1), it yields
(2)t′1=t1,(3)t′m+1=inttm+1−t1∗p−1∗p+t1,where t1′ and tm+1′ are the beginning time and the ending time after interpolation, respectively; the value of inttm+1−t1∗p−1 is an integer multiple relationship between the difference of time tm+1−t1 and the interpolation interval p. So the time range after interpolation is t1′,tm+1′.Similarly, set the same interpolation interval to interpolate the second segment data fromtm+1,ym+1 to t2m+1,y2m+1 in Y. The time range of the original data is tm+1,t2m+1.From formula (1), we can obtain
(4)t′m+2=t′m+1+p,(5)t′2m+1=intt2m+1−t′m+2∗p−1∗p+t′m+2,where tm+2′tm+2′>tm+1 and t2m+1′t2m+1′>t2m+1 are the beginning time and the ending time after interpolation, respectively.According to (2), (3), (4) and (5), we can obtain
(6)t′i+1=t′i+p,(7)t′m+i=inttm+i−t′i+1∗p−1∗p+t′i+1.Based on the above analysis, the algorithm above is described as follows:Step 1.
Wheni=1, the observed data is y=t1,y1,t2,y2,⋯,tm,ym,tm+1,ym+1. The beginning time of interpolation data and the observed data are the same.Step 2.
The interpolation time data based on formulas (2) and (3), and it is equal.Step 3.
Lagrange Interpolation method is used to interpolate the MW speed data by formula (1).Step 4.
Wheni>1, the observed data is y=ti,yi,ti+1,yi+1,⋯,ti−1+m,yi−1+m,ti+m,yi+m.Step 5.
The interpolation time data based on formulas (6) and (7), and it is also equal.Step 6.
The Lagrange interpolation method is also used to interpolate the MW speed data by formula (1).Step 7.
Ifi≠N−m, back to Step 4. Finally, the equal interval data is expressed by Y′=t1′,y1′,t2′,y2′,⋯,tN1′,yN1′T.The flow chart of the equal interval interpolation method is shown in Figure2.Figure 2
The flow chart of equal interval interpolation method.


### 2.3. Determination of the Period of the Momentum Wheel Speed Data

The equal interval data can be calculated by the Lagrange interpolation method, and it is approximately periodic. In order to determine the period of the MW speed data, the Fast Fourier Transformation (FFT) method is adopted to calculate the frequency of the equal interval dataY′=t1′,y1′,t2′,y2′,⋯,tN′,yN′T. If the calculated frequency is a constant, it shows that the MW speed data has the characteristic of approximate periodicity.Setting the time interva1p=0.1s, then the sampling frequency is fs=1/p=10Hz. The FFT method [10] is used to analyze the frequency of the MW speed data. The frequency diagram of MW data is shown in Figure 3.Figure 3
The frequency analysis of the momentum wheel speed.According to the analysis of data and Figure3, the amplitude of the 35th data is maximum, so the frequency of MW speed data f=0.0416Hz and the period of the MW speed T=f−1=24.02s≈24s.


## 2.1. Analysis of the Momentum Wheel Speed Data

The momentum wheel speed data from a simulation test of communication satellite is time series data. Assuming that the observed MW speed isyti at time ti, written as yi, the change curve of MW speed is shown in Figure 1. Because of accidental factors such as cosmic ray and accidental contact failure of circuit, the sensor data is abnormal [7]; for example, the observed data is incompletely equal interval distribution and the partially observed data loss. Meanwhile, it can be seen that the MW speed is approximately periodic and its amplitude grows linearly. In order to predict MW speed, the Lagrange interpolation method [8] is used to calculate the MW speed, so that the unequal interval data is converted into equal interval data.Figure 1
The change curves of momentum wheel speed.


## 2.2. Preprocessing of Equal Interval Interpolation of Momentum Wheel Speed Data

The time series data of the observed MW speed isY=t1,y1,t2,y2,⋯,tN,yNT, where N is the total number of data. However, Y is the unequal interval time series data, and the forecast model based on equal interval time series data. So Y has to convert into equal interval time series data by the interpolation method.Firstly, set the interpolation intervalp. To interpolate the data segment from t1,y1 to tm+1,ym+1 in Y, the Lagrange interpolation formula [9] can be given by
(1)Lmt=∑k=1m+1∏i=0i≠0m+1t−titk−tiyk,where t is the interpolation time; t1,t2,⋯,tm+1 are the observation time which satisfy t1<t2<⋯<tm+1.According to formula (1), it yields
(2)t′1=t1,(3)t′m+1=inttm+1−t1∗p−1∗p+t1,where t1′ and tm+1′ are the beginning time and the ending time after interpolation, respectively; the value of inttm+1−t1∗p−1 is an integer multiple relationship between the difference of time tm+1−t1 and the interpolation interval p. So the time range after interpolation is t1′,tm+1′.Similarly, set the same interpolation interval to interpolate the second segment data fromtm+1,ym+1 to t2m+1,y2m+1 in Y. The time range of the original data is tm+1,t2m+1.From formula (1), we can obtain
(4)t′m+2=t′m+1+p,(5)t′2m+1=intt2m+1−t′m+2∗p−1∗p+t′m+2,where tm+2′tm+2′>tm+1 and t2m+1′t2m+1′>t2m+1 are the beginning time and the ending time after interpolation, respectively.According to (2), (3), (4) and (5), we can obtain
(6)t′i+1=t′i+p,(7)t′m+i=inttm+i−t′i+1∗p−1∗p+t′i+1.Based on the above analysis, the algorithm above is described as follows:Step 1.
Wheni=1, the observed data is y=t1,y1,t2,y2,⋯,tm,ym,tm+1,ym+1. The beginning time of interpolation data and the observed data are the same.Step 2.
The interpolation time data based on formulas (2) and (3), and it is equal.Step 3.
Lagrange Interpolation method is used to interpolate the MW speed data by formula (1).Step 4.
Wheni>1, the observed data is y=ti,yi,ti+1,yi+1,⋯,ti−1+m,yi−1+m,ti+m,yi+m.Step 5.
The interpolation time data based on formulas (6) and (7), and it is also equal.Step 6.
The Lagrange interpolation method is also used to interpolate the MW speed data by formula (1).Step 7.
Ifi≠N−m, back to Step 4. Finally, the equal interval data is expressed by Y′=t1′,y1′,t2′,y2′,⋯,tN1′,yN1′T.The flow chart of the equal interval interpolation method is shown in Figure2.Figure 2
The flow chart of equal interval interpolation method.


## 2.3. Determination of the Period of the Momentum Wheel Speed Data

The equal interval data can be calculated by the Lagrange interpolation method, and it is approximately periodic. In order to determine the period of the MW speed data, the Fast Fourier Transformation (FFT) method is adopted to calculate the frequency of the equal interval dataY′=t1′,y1′,t2′,y2′,⋯,tN′,yN′T. If the calculated frequency is a constant, it shows that the MW speed data has the characteristic of approximate periodicity.Setting the time interva1p=0.1s, then the sampling frequency is fs=1/p=10Hz. The FFT method [10] is used to analyze the frequency of the MW speed data. The frequency diagram of MW data is shown in Figure 3.Figure 3
The frequency analysis of the momentum wheel speed.According to the analysis of data and Figure3, the amplitude of the 35th data is maximum, so the frequency of MW speed data f=0.0416Hz and the period of the MW speed T=f−1=24.02s≈24s.


## 3. Modeling and Prediction of Momentum Wheel Data

### 3.1. Determination of the Order of the MW Speed Data

The ARMA model is established for the equal interval dataY′=t1′,y1′,t2′,y2′,⋯,tN′,yN′T, and the order of the model is determined by general AIC (Akaike Information Criterion) and BIC (Bayes Information Criterion). The loss functions of AIC and BIC [11] are shown as follows:
(8)EAIC=lnδ2+2i+jn,EBIC=lnδ2+2nn−i−1i+jn,where σ2 is the variance of error between the model and the observed data, i and j are the order of the AR and MA model, respectively, and n is the length of data sequence for modeling and satisfies n=5∗i+j+1.The loss value of AIC and BIC criterion under the lengthn of the data sequence is shown in Figure 4. We can see that the loss function values of AIC and BIC are relatively smoother between n=10 and n=18, which means the robustness of the ARMA model is better. When n=16, the loss values of AIC and BIC are minimum. Combining n=5∗i+j+1 with i>j, it can be calculated when i=3 and j=0, and the following study will be based on the third order AR model.Figure 4
The loss value of AIC and BIC criterion under the lengthn of the data sequence.


### 3.2. AR Modeling of the Momentum Wheel Speed Data

In order to extract the MW speed data conveniently, we convert the equal interval data sequenceY′=t1′,y1′,t2′,y2′,⋯,tN′,yN′T into Y″=1,y1′,2,y2′,⋯,N,yN′T, and the period of Y″ is calculated as T′=T/p=240. The curve of MW speed data is similar with the sinusoidal function. In order to reduce the calculation, the data of 0, quarter, a half, and three quarter periods are selected as characteristic points expressed by h1, h60, h120, and h180. We find that each of them is basically linear in Figure 5. Therefore, an AR model can be established for each element [12–14].Figure 5
The curves of the selected data.As discussed above, the loss function is minimized with a third order AR model, the AR equation is shown as follows:(9)xk=a1xk−1+a2xk−2+a3xk−3+ek,where a1, a2, and a3 are the AR model coefficients, ek is a white Gaussian noise with zero mean and variance δe2=0.09, xkk=0,1,2,⋯,intN1/T′ is hi=yi′,yi+T′′,yi+2T′′,⋯,yi+k∗T′′, respectively.AssumeφTk=xk−1,xk−2,xk−3 and θ=a1,a2,a3T, the formula (9) can be written as a vector form by
(10)xk=φTkθ+ek.According to (10), we can obtain
(11)ek=xk−φTkθ.The least-square method is used to calculate the minimum of RSS (residual sum of squares) for unknown parameter vectorθ. The performance function Jθ is
(12)Jθ=∑i=1te2k=∑i=1kxi−φTiθ2.The purpose of equation (12) is to obtain the estimation value θ^, when the function Jθ is minimized.Formula (10) is rewritten as a vector form
(13)zk=Hkθ+vk,where
(14)zk=x1x2⋮xk,Hk=φT1φT2⋮φTk,vk=e1e2⋮ek.The least-square estimation ofθ [15] at time k is
(15)θ^k=HTkHk−1HTkzk.


### 3.3. Verification of the MW Speed AR Model

The first 6000 data is selected as a training dataset, and others is as a test data; the estimated value of model coefficient can be obtained by equation (15).
(16)θ^k=a^1,a^2,a^3T.Based on formulas (10) and (15), the value of xk is obtained which can be used to predict the value of xk+1. Therefore, the MW speed data in the next period can be obtained by the AR model. The forecast model is shown as follows:
(17)x^k+1=φTk+1θ^k+ek+1.The predicted value and the measured value are shown in Figure6. It can be seen that the predicted data is in accordance with the measured data. The error between the measured value and the predicted value is shown in Figure 7. The results show the effectiveness and feasibility of the proposed algorithm.Figure 6
The prediction value and measured value.Figure 7
The error between the predicted value and measured value.The standard deviation of the residuals is used to verify the accuracy of the forecast model, and it is expressed by(18)δ=∑i=1T′xi−x^i2T′,i=1,2,⋯,T′,where xi is the actual observed data, and x^i is the predicted data of the forecast model. The results are shown in Table 1.Table 1
The standard deviation of the residual between the measured data and the predicted data.
Order
First order
Second order
Third order
Fourth order
The standard deviation of the residualδ
0.4052
0.0599
0.0575
0.1214From Table1, it can be seen that the standard deviation of the residual is minimized when the order of the forecast model is three. So the third order forecast model is more accurate than others.


## 3.1. Determination of the Order of the MW Speed Data

The ARMA model is established for the equal interval dataY′=t1′,y1′,t2′,y2′,⋯,tN′,yN′T, and the order of the model is determined by general AIC (Akaike Information Criterion) and BIC (Bayes Information Criterion). The loss functions of AIC and BIC [11] are shown as follows:
(8)EAIC=lnδ2+2i+jn,EBIC=lnδ2+2nn−i−1i+jn,where σ2 is the variance of error between the model and the observed data, i and j are the order of the AR and MA model, respectively, and n is the length of data sequence for modeling and satisfies n=5∗i+j+1.The loss value of AIC and BIC criterion under the lengthn of the data sequence is shown in Figure 4. We can see that the loss function values of AIC and BIC are relatively smoother between n=10 and n=18, which means the robustness of the ARMA model is better. When n=16, the loss values of AIC and BIC are minimum. Combining n=5∗i+j+1 with i>j, it can be calculated when i=3 and j=0, and the following study will be based on the third order AR model.Figure 4
The loss value of AIC and BIC criterion under the lengthn of the data sequence.


## 3.2. AR Modeling of the Momentum Wheel Speed Data

In order to extract the MW speed data conveniently, we convert the equal interval data sequenceY′=t1′,y1′,t2′,y2′,⋯,tN′,yN′T into Y″=1,y1′,2,y2′,⋯,N,yN′T, and the period of Y″ is calculated as T′=T/p=240. The curve of MW speed data is similar with the sinusoidal function. In order to reduce the calculation, the data of 0, quarter, a half, and three quarter periods are selected as characteristic points expressed by h1, h60, h120, and h180. We find that each of them is basically linear in Figure 5. Therefore, an AR model can be established for each element [12–14].Figure 5
The curves of the selected data.As discussed above, the loss function is minimized with a third order AR model, the AR equation is shown as follows:(9)xk=a1xk−1+a2xk−2+a3xk−3+ek,where a1, a2, and a3 are the AR model coefficients, ek is a white Gaussian noise with zero mean and variance δe2=0.09, xkk=0,1,2,⋯,intN1/T′ is hi=yi′,yi+T′′,yi+2T′′,⋯,yi+k∗T′′, respectively.AssumeφTk=xk−1,xk−2,xk−3 and θ=a1,a2,a3T, the formula (9) can be written as a vector form by
(10)xk=φTkθ+ek.According to (10), we can obtain
(11)ek=xk−φTkθ.The least-square method is used to calculate the minimum of RSS (residual sum of squares) for unknown parameter vectorθ. The performance function Jθ is
(12)Jθ=∑i=1te2k=∑i=1kxi−φTiθ2.The purpose of equation (12) is to obtain the estimation value θ^, when the function Jθ is minimized.Formula (10) is rewritten as a vector form
(13)zk=Hkθ+vk,where
(14)zk=x1x2⋮xk,Hk=φT1φT2⋮φTk,vk=e1e2⋮ek.The least-square estimation ofθ [15] at time k is
(15)θ^k=HTkHk−1HTkzk.


## 3.3. Verification of the MW Speed AR Model

The first 6000 data is selected as a training dataset, and others is as a test data; the estimated value of model coefficient can be obtained by equation (15).
(16)θ^k=a^1,a^2,a^3T.Based on formulas (10) and (15), the value of xk is obtained which can be used to predict the value of xk+1. Therefore, the MW speed data in the next period can be obtained by the AR model. The forecast model is shown as follows:
(17)x^k+1=φTk+1θ^k+ek+1.The predicted value and the measured value are shown in Figure6. It can be seen that the predicted data is in accordance with the measured data. The error between the measured value and the predicted value is shown in Figure 7. The results show the effectiveness and feasibility of the proposed algorithm.Figure 6
The prediction value and measured value.Figure 7
The error between the predicted value and measured value.The standard deviation of the residuals is used to verify the accuracy of the forecast model, and it is expressed by(18)δ=∑i=1T′xi−x^i2T′,i=1,2,⋯,T′,where xi is the actual observed data, and x^i is the predicted data of the forecast model. The results are shown in Table 1.Table 1
The standard deviation of the residual between the measured data and the predicted data.
Order
First order
Second order
Third order
Fourth order
The standard deviation of the residualδ
0.4052
0.0599
0.0575
0.1214From Table1, it can be seen that the standard deviation of the residual is minimized when the order of the forecast model is three. So the third order forecast model is more accurate than others.


## 4. Conclusion

This paper presents the MW speed AR model based on the variation law of MW speed of the three-axis stabilized GEO satellite on the orbit. It can be used to estimate the future data trends. The model will play an important role in the actual control of satellite.By analyzing the data variation of MW rotation rate, the Lagrange interpolation method is proposed to solve the problem of unequal interval and data loss. Analyzing the data feature of the MW, the FFT algorithm is used to calculate the period of the MW speed. Establish the multidimensional AR model to verify the accuracy of the model. The model order and prediction results can be obtained.The accuracy of the model prediction is compared by variances. The results show that the third order AR multidimensional model has a higher performance.



---

*Source: 5142696-2020-03-03.xml*