A novel chaotic time series prediction method and its application to carrier vibration interference attitude prediction of stabilized platform

2.1. GM(1, 1) model

Assume $X^{\left(0\right)}$ is a nonnegative time sequence:

$X^{\left(1\right)}$ is the 1-AGO (1-accumulated generating operation) sequence of $X^{\left(0\right)}$:

where, $x^{\left(1\right)}\left(k\right)={\sum }_{i=1}^k{x}^{\left(0\right)}\left(i\right)$, $k=\text{1,}\text{2,...}\text{,}n$.

The sequence mean generation of consecutive neighbors of $X^{\left(1\right)}$ is $Z^{\left(1\right)}$:

where, $z^{\left(1\right)}\left(k\right)=\left(x^{\left(1\right)}\right(k)+x^{\left(1\right)}(k-1\left)\right)/2$, $k=\text{2,}\text{3,...}\text{,}n$.

Definition 1. Let $x^{\left(0\right)}\left(k\right)$, $x^{\left(1\right)}\left(k\right)$ and $z^{\left(1\right)}\left(k\right)$ be as above mentioned, and the grey differential equation:

is called GM(1, 1) Model.

Theorem 1. Based on Eqs. (1)-(3), if $\widehat{a}=\left(\begin{array}{l}a\\ b\end{array}\right)$ is the model parameter, and:

So the model parameters evaluated by least square method of grey differential equation $x^{\left(0\right)}\left(k\right)+a{z}^{\left(1\right)}\left(k\right)=b$ satisfy:

Let$x^{\left(0\right)}\left(k\right)\text{,}$$x^{\left(1\right)}\left(k\right)$ and $z^{\left(1\right)}\left(k\right)$ be as above mentioned, Eq. (7) is called whitening differential equation of GM(1, 1):

Assume $B$, $Y$ and $\widehat{a}$ is described in Theorem 1 and $\widehat{a}=\left(\begin{array}{l}a\\ b\end{array}\right)=(B^{T}B{)}^{-1}{B}^{T}Y$, and apply the inverse AGO (IAGO) to figure out ${\widehat{x}}^{\left(0\right)}(k+1)$, so we can obtain:

And ${\widehat{x}}^{\left(0\right)}\left(1\right)={\widehat{x}}^{\left(1\right)}\left(1\right)=x^{\left(0\right)}\left(1\right)$.

2.2. Modifications of GM(1, 1) (MGM(1, 1))

As is known to all, many improvements on GM(1, 1) model have been published and achieved. What’s more, there are three main approaches: optimizing background value, optimizing original condition and modifying parameter estimation method. On the basis of GM(1, 1), we draw lessons from modifications of GM(1, 1) model to optimize original condition and background value of GM(1, 1) meanwhile in this paper.

3.1. Phase space reconstruction technique

In order to reveal hidden patterns and trends in chaotic time series, phase space reconstruction technique is always employed. And currently the most widely used approach is phase space delay coordinate reconstitution theory. Phase space reconstruction is also called dynamic system reconstruction, and the most popular approach is phase space delay coordinate reconstitution theory which holds that any component evolution in the system is determined by other correlative components, and the most important characteristics of original phase space portrait can be kept in the state trajectory derived by embedding single variable time series to a new coordinate system. Let the one-dimension time series be denoted as $\mathbf{X}=\left\{x\right(k),k=1,\mathrm{}2,\mathrm{}\dots ,\mathrm{}n\}$, and apply phase space delay coordinate reconstitution method to construct a high-dimension space $\mathbf{Y}$ with a given time delay $\tau$ and a certain embedding dimension $m$:

where $M=n-(m-1)\tau$ is the number of points in the phase space.

3.2. C-C method

Selections of time delay $\tau$ and embedding dimension $m$ are very important in phase space reconstruction. At present, there are two main viewpoints about selections of time delay $\tau$ and embedding dimension $m$. One is that the two parameters should be selected independently because they are non-correlative. The other one holds the opposite view. C-C method [28] witch believes that $m$ and $\tau$ are correlative is applied to calculate the two parameters with the help of correlation integral in this paper. The following function is the correlation integral for the embedded time series:

where $\mathrm{\Theta }\left(a\right)=0$, if $a<0$, $\mathrm{\Theta }\left(a\right)=1$, if $a\ge 0$, and the correlation dimension is defined as:

Brock [29] studied the BDS statistic, which is based on the correlation integral, to test the null hypothesis that a given data set is independently and identically distributed (iid). The BDS statistic originates from the statistical properties of the correlation integral, and it measures the statistical significance of calculations of the correlation dimension. The BDS statistic defined by:

With the help of $C(1,r)$ and ${\sigma }^2(m,r)$ calculated respectively by $C(1,n,r,t)$ and ${\widehat{\sigma }}^2$, the BDS statistic becomes:

And this converges to a standard normal distribution as $n\to \mathrm{\infty }$.

Where:

In order to study the nonlinear dependence and eliminate spurious temporal correlations, we must subdivide the time series $\mathbf{X}$ into $t$ disjoint time series. For general $t$, $S(m,r,t)$ can be further computed by:

Finally, as $n\to \mathrm{\infty }$, we can write:

We select several representative values $r_j$, and define the quantity:

which is a measure of the variation of $S(m,r,t)$ with $r$. The locally optimal times $t$ are then the zero crossings of $S(m,r,t)$ and the minima of $\mathrm{\Delta }S(m,t)$. The zero crossings of $S(m,r,t)$ should be nearly the same for all $m$ and $r$, and the minima of $\mathrm{\Delta }S(m,t)$ should be nearly the same for all $m$ (otherwise, the time is not locally optimal). The delay time $\tau$ will correspond to the first of these locally optimal times.

3.3. Identification of chaos time series

Noise is unavoidable in measured data. But different from random noises, Chaos is the pseudorandom behavior that appears in confirmed non-linear systems, and should not be confused with random noises. Therefore, before analyzing time series with chaotic methods, identification of chaos should be finished first. Now there have been many methods to identify chaos, such as the Power spectrum method, the method of principal component analysis, the method of Poincare surface of section, the Lyapunov exponent method and so on.

Power spectrum method [30] is mainly to analyze the power spectrum of time series. If the spectrum has a single peak (or several peaks), it corresponds to the periodic (or quasi periodic) sequence; if there is no obvious peak or peak to be connected, then corresponds to chaotic sequence. Periodic motion in the power spectrum corresponds to peak, and the characteristic of chaos is the appearance of noise background and wide peaks in the spectrum. But the method has a large amount of calculation and the performance requirements are also high too. Principal component analysis (PCA) [30] is a new method to identify chaos and noise in recent years. Because there is a big difference between the noise and the chaotic signal in the main components, the main component spectrum of the noise is a straight line which is parallel to the $X$ axis, and the main component spectrum of the chaotic signal is a line witch pass a fixed point and the slope of is negative. Therefore, the standard deviation of principal component distribution can be used as a feature of distinguishing chaos and noise. This method needs to reconstruct the phase space and carry out a large number of mathematical matrix, then analyze the shape of the main component spectrum, and the procedure is more complicated. Poincare section method [30] is suitable to select a section in the phase space, and this section is usually referred to as the Poincare cross section. The cross section of continuous trajectory of the phase space and the intersection of the Poincare is called the cut-off point. When the Poincare section is a dense point with a split type structure, it is chaotic; when the Poincare cross section has only one or a few discrete points, or a closed curve, it is periodic motion or quasi periodic motion. Although the method is simple to determine, but for Poincare section selection has a certain degree of difficulty, such as not with the section line tangent, more cannot contain trajectory and so on.

Lyapunov exponent which is the long-term average result along the trajectory is a global feature. The Lyapunov exponent is an important parameter for depicting the characteristic of nonlinear time series, and it is the most basic and useful dynamical parameter for deterministic chaotic system. If the largest Lyapunov exponent is less than zero, the system evolution must have a confirmed rule, and it isn’t a chaotic system, vice versa. In the paper, the Lyapunov exponent method will be employed to identify whether the measured series from stabilized platform experiment is a chaotic series or not.

The Wolf method can be used to compute the largest Lyapunov exponent [31]. Take $Y\left(t_0\right)$ as the initial point, and the distance between $Y\left(t_0\right)$ and its nearest point $Y_0\left(t_0\right)$ is $L_0$. Measure the divergence of the two points, if their distance meets $L_0^{\text{'}}=\left|Y\left(t_1\right)-Y_0\left(t_1\right)\right|>\epsilon$$(\epsilon >0)$ at time $t_1$, find another point $Y_1\left(t_1\right)$ around $Y\left(t_1\right)$ which satisfies $L_1=\left|Y\left(t_1\right)-Y_1\left(t_1\right)\right|<\epsilon$ and the angle between $Y_1\left(t_1\right)$ and $Y\left(t_1\right)$ is small enough. Continue the above process until $Y\left(t\right)$ reaches the end of the sequence. In the process, the total iteration times is $M$, and the largest Lyapunov exponent can be defined as:

3.4. Chaos prediction with phase point L1 norm

In above paragraphs, we have constructed a high-dimension phase space $\mathbf{Y}$ for raw time series $\mathbf{X}=\left\{x\right(k),\mathrm{}k=1,\mathrm{}2,\dots ,\mathrm{}n\}$ with time delay $\tau$ and embedding dimension $m$. Define $L_1$ norm of phase point $\mathbf{Y}\left(k\right)$ at time $k$ as [32]:

and $L_1$ norm series ${\mathbf{L}}^{\left(0\right)}=\left\{l^{\left(0\right)}\right(k),\mathrm{}k=1,\mathrm{}2,\dots ,\mathrm{}N\}$of phase points can be obtained, where $N=n-(m-1)\tau$. From Eq. (26) we can find that, $L_1$ norm series ${\mathbf{L}}^{\left(0\right)}$ is a non-negative series. And it needs to ensure the class ratio $h^{\left(0\right)}\left(k\right)$ of ${\mathbf{L}}^{\left(0\right)}$ to be within the interval $(e^{-2/(N+1)},\mathrm{}{e}^{2/(N+1)})$ [33]:

Then the MGM(1, 1) model presented in the next section could be used to make predictions for $L_1$ norm series, With the predicted value ${\widehat{l}}^{\left(0\right)}\left(N+1\right)\text{,}$ the following equation can be established according to L1 norm definition:

In Eq. (28), only $x(N+1+(m-1\left)\tau \right)$ is unknown, and $x(N+1+(m-1\left)\tau \right)$ is just the value to be predicted. Therefore, solving Eq. (28) can obtain the predicted value $\widehat{x}(n+1)$ of raw time series:

5.1. The stabilized platform and the sample data

For showing the prediction capability and efficiency of the proposed MGM(1, 1) prediction model, take an actual stabilized platform system as an example to analyze. The stabilized platform is shown in Fig. 2. The stabilized platform consists of two level platforms witch include the bottom stabilized platform and the upper tracking platform. The bottom stabilization platform is mainly to reduce the disturbance of the external disturbance of the base, and provides a good level reference for target tracking. After ensuring the level of the base, the upper tracking platform can realize the effective tracking of the target. The stable platform system is drived by servo motor. Firstly, the adjusting shafts (pitch and yaw axis) of the bottom stabilized platform are contained by the main control computer with the control information (after calculation and processing of attitude adjustment information) of regulation to ensure stable layer approximation level. Secondly, the upper tracking platform according to the information of goal drives the servo motor for rotation of the azimuth and pitch regulation so as to achieve the goal of fast tracking. What’s more, the stable platform load is $m\approx \mathrm{}$30 kg, the bottom stabilized platform diameter is $d=\mathrm{}$450 mm, the frequency of stable platform is $f=\mathrm{}$1 Hz, the stable layer regulating shaft adjusting range is ±20°. The range of azimuth adjustment and pitch angle range of the upper tracking platform is 360°and ±5°respectively.

In Fig. 3, Signals ${\theta }_i$, ${\theta }_o$, $e$ denote the reference input, system output and tracking error of stabilized platform control system respectively. $G_{APR}\left(s\right)$ and $G_{rob}\left(s\right)$ is the regulator and object of position control respectively. $e^{-\tau s}$ is the link of time, ${\theta }_d$ is the carrier interference attitude, $D\left(s\right)$ is the interference gallery transfer function and $CPM$ is the chaos time prediction module. The stabilized platform will be influenced by the carrier vibration interference attitude. Therefore, if the information of carrier vibration interference attitude can be obtained by prediction, the stabilized platform control system can suppress vibration interference in advance based on interference suppression. For this reason, the carrier vibration interference attitude prediction is very significant.

Fig. 2The stabilized platform

Fig. 3The compound control system structure

The sample data of carrier interference attitude are measured at a sampling frequency of 50 Hz with an Innalabs Attitude and Heading Reference System (AHRS) in a stabilized platform experiment. The data length is 1000. Take pitch attitude to analyze, the former 850 data are employed as training samples and latter 150 data are used as testing samples. The raw sample time series is shown in Fig. 4, the former 850 data and the latter 150 data are shown in Fig. 5 and Fig. 6.

Fig. 4Raw time series of pitch attitude

Fig. 5The former 850 points of sample data

Fig. 6The latter 150 points of sample data

Before reconstructing the phase space, the pitch attitude series should be judged by the largest Lyapunov exponent method. Table 1 shows the largest Lyapunov exponent $\lambda$ is greater than zero to prove that the pitch attitude series is chaotic and can be analyzed with chaotic methods.

By reconstructing the phase space $R^d$ with time delay $\mathrm{\tau }$ and embedding dimension $m$, the L1 norm series of phase points can be obtained from Eq. (28) and the curve of L1 norm series is illustrated in Fig. 7. From Figs. 4 and 7, we can see that raw time series and L1 norm series are both oscillatory. Before predicting, class ratios should be examined. Since the raw time series of pitch attitude is not a nonnegative time sequence, it is preprocessed with upward translation in order to make the L1 norm series within the interval $(e^{-2/(N+1)},\mathrm{}{e}^{2/(N+1)})$ [33] as is shown in Fig. 6. The latest five points are employed to construct grey models, and the interval of class ratios would be $(e^{-2/(5+1)},\mathrm{}{e}^{2/(5+1)})=\text{(0.7165,}\text{1.3956)}$. From Fig. 8 we can see that, the class ratios of roll attitude series after upward translation and L1 norm series all lie in the interval $(e^{-2/(5+1)},\mathrm{}\mathrm{}{e}^{2/(5+1)})=\text{(0.7165,}\text{1.3956)}$. In order to compare the performances of GM(1, 1) model, residual GM(1, 1) mode (RGM) [6] and the modified GM(1, 1) model (MGM) for oscillatory sequences, firstly of all, the three models are formed with raw time series respectively. Then five other chaotic prediction methods, Tang’s method (method 1), adding-weight one-rank local-region method (method 2), GM(1, 1) model applied in reconstructed phase space (method 3), residual GM(1, 1) model applied in reconstructed phase space(method 4), the combined forecasting model (method 5) [26] and the proposed MGM(1, 1) model applied in reconstructed phase space (method 6) proposed in the paper are used to predict the latter 150 sample points respectively.

Fig. 7L1 norm series

Fig. 8Class ratios of L1 norm series

5.2. Simulation results

In comparisons with the actual pitch data series curves, Figs. 9, 10 and 11 depict the predicted pitch attitude curves using GM(1, 1) model, residual GM(1, 1) model and MGM(1, 1) model respectively. And in order to compare the prediction method proposed in the paper with other mentioned methods, three criteria that contain the mean absolute percent error (MAPE), root mean square error (RMSE) and forecasting effective measure(FEM) are employed here [6].

Form Fig. 9 we can see that, big errors always exist at the inflexions by predicting with GM(1, 1) model. In Fig. 10, the large errors have been controlled when predicting with the residual GM(1, 1) model in a certain extent. On the contrary, when predicting with MGM in Fig. 11, errors can be reduced to a large certain extent at the inflexions due to optimizing background value and original condition. The prediction evaluation index is shown in Table 2, the overall performance of MGM is superior to GM(1, 1) model and residual GM(1, 1) model. Therefore we can conclude that, the MGM(1, 1) model would have a better performance for oscillatory sequences.

Fig. 9Prediction curve of GM(1, 1)

Fig. 10Prediction curve of RGM(1, 1)

Fig. 11Prediction curve of MGM(1, 1)

Figs. 12-17 are prediction curves of each model, and Table 2 lists the evaluation results.

Figs. 12-17 depict predictions made from reconstructed phase space. From above figures, we can see that method 1 has the worst prediction performance in Fig. 12, and the evaluation results indicate that it is not fit for the attitude prediction of stabilized platform. Fig. 13 indicates that the performance of adding-weight one-rank local-region method is quite unstable because it is a simple prediction method. Compare Fig. 12 with Fig. 13, we can find that method 2 is better than method 1. AOLM finds a set of neighbor points in reconstructed phase space, and constructs a first-order linear model. Then adopt least square method to estimate the model coefficients with these neighbor points, and system prediction model would be obtained in this way. The poor capability of AOLM may be due to its first-order linear model, hence choosing a proper model in reconstructed phase space is also very important. When comparing method 6 with method 1 and method 2, we can see that there is no bad prediction inflexion in Fig. 16 which exists in Fig. 12 and Fig. 13. And Table 2 shows that, MAPE and RMSE of method 5 are almost decuple better than the method 1 and method 2, and FEM is high to 0.97986.

Fig. 12Prediction curve of Tang’s method (method 1)

Fig. 13Prediction curve of AOLM (method 2)

Fig. 14Prediction curve of GM-RPS (method 3)

Fig. 15Prediction curve of RGM-RPS (method 4)

Fig. 16Prediction curve of CFM (method 5)

Fig. 17Prediction curve of MGM-RPS (method 6)

By comparing Figs. 14, 15 and 16 with Figs. 9, 10 and 11, we can find, that predictions made in reconstructed phase space substantially reduce the errors of GM(1, 1), residual GM(1, 1) model and MGM predictors. GM-RPS is better than GM(1, 1) obviously, and evaluation results show that its performance is almost equivalent to MGM. What’s more, RGM-RPS is superior to RGM and GM distinctly. By comparing Fig. 14, with Fig. 15, the RGM-RPS is better than GM-RPS, because the RGM-RPS has compensated the predicted value by counting residual series ${\mathbf{e}}^{\left(0\right)}$ to improve the prediction precision. By comparing Figs. 14, 15 with Fig. 16 respectively, among the three methods, MGM-RPS is the best with no doubt. MGM-RPS has not only compensated the predicted value, but also optimized background value and optimized original condition. Furthermore, in Table 2, we can see that MAPE and RMSE of method 5 are almost decuple better than the method 3 and method 4 too. MAPE and RMSE of method 3 are 0.020832 and 0.1634, and MAPE and RMSE of method 4 are 0.017854 and 0.1013 respectively. However, the MAPE and RMSE of method 5 are only 0.014764 and 0.0938.

When comparing Fig. 16 with Fig. 17, the prediction performances shown in Table 2 exhibit that, the MAPE and RMSE of method 5 is bigger than the method 6. And the computational cost of method 5 is high. The overall performance of MGM-RPS is distinctly superior to combined forecasting model. Consequently we can reach a conclusion with no doubt: the newly proposed chaos prediction method with MGM(1, 1) model applied in reconstructed phase space could effectively extract the chaotic time series features.