CN105320845A - Time sequence forecast method based on quantum gravity algorithm - Google Patents

Abstract

The invention discloses a time sequence forecast method based on a quantum gravity algorithm. A time sequence forecast model is constructed by adopting a Volterra series expanded formula. The method comprises the steps that firstly, an input signal vector subjected to Volterra series P-order truncation is constructed on the basis of quoting a phase space framing reconstruction technology; secondly, acceleration is introduced to serve as a variable parameter into the quantum gravity algorithm put forward by Mohadeseh Soleimanpour et al., and a Volterra kernel function of the forecast model is effectively trained by means of the algorithm; finally, a forecast value of forecast point time is obtained through a linear combination of the input signal vector of the forecast point time and the Volterra kernel function. It is verified through experiments that the accuracy of time sequence forecast can be effectively improved by introducing the quantum gravity algorithm into identification of the Volterra kernel function of a non-linear system.

Description

A kind of Time Series Forecasting Methods based on quantum gravity algorithm

Technical field

The invention belongs to fields such as signal transacting and signal estimations, more specifically, it is related to a kind of Time Series Forecasting Methods based on quantum gravity algorithm.

Background technology

Forecast of Nonlinear Time Series method develops since the 1980s mid-term, obtained in recent years deeper into research and widely application.Wherein, Volterra series models analysis method because its output be its filter kernel linear function, its filtering performance can be analyzed with existing linear tool, it has also become one of widely used forecast model.Volterra series is a kind of functional, and most of nonlinear dynamic systems can use Volterra series approximations to any accurately degree.In other words, the nonlinear dynamic system of a major class has its intrinsic Volterra kernel function.Volterra kernel functions can characterize the property of nonlinear system completely.Therefore, only it is to be understood that core is just enough to determine response of the system to any input stimulus.But, nonlinear system Volterra series models are set up, are not an easy things, the kernel function estimation of its high-order is the maximum difficulty that it faces, and this largely constrains effective application of Volterra functional series models.And existing Volterra kernel functions identification mainly uses traditional least square method, the method is scanned for using gradient information, is easily absorbed in Local Minimum, cannot get satisfactory result, especially in the case of serious interference, and the stability of identification is unsatisfactory.Its result has led to the inaccurate of modeling, influence the estimated performance of model, especially training effect is easily influenceed when the information related to predicted orbit included in training set is very few or is mingled with much unrelated information, cause model training insufficient or over-fitting, cause multi-step prediction performance not good.

And quantum gravity algorithm (QuantumInspiredGravitationalSearchAlgorithm, abbreviation QGSA) it is used as a kind of new intelligent optimization algorithm, its improved gravitation searching algorithm (GravitationalSearchAlgorithm proposed from quantum mechanics angle, abbreviation GSA) model, particle is defined in a vector subspace, the vector subspace is determined (particle in Two-dimensional hole gas can appear in any spatial point with certain probability) by probability density function, think that particle has the behavior of quantum, it can be scanned in whole solution space, particular problem to demand solution is not required particularly in itself, thus ability of searching optimum is stronger, overcome some shortcomings of above-mentioned Volterra kernel functions discrimination method, in addition, also there is easily realization, the advantages of convergence is soon and parameter is few.Therefore, Volterra kernel function identification and its Time Series Forecasting Methods of the research based on quantum gravity algorithm has important theory value and realistic meaning.

The content of the invention

It is an object of the invention to overcome the deficiencies in the prior art, a kind of Time Series Forecasting Methods based on quantum gravity algorithm are provided, quantum gravity algorithm is introduced into the identification of the Volterra kernel functions of nonlinear system, to improve the degree of accuracy of electronic product degenerate state prediction.

For achieving the above object, the Time Series Forecasting Methods of the invention based on quantum gravity algorithm, it is characterised in that comprise the following steps：

(1), data prediction

(1.1) the data set DS={ x (1), x (2) ..., x (D) } of pre-selection, is set, wherein D is data amount check；

(1.2) data set DS smallest embedding dimension number d and delay time T, is determined according to the method for differential entropy rate；

(1.3) framing reconstruct, is carried out to data set DS using the function windowize in Matlab, and each frame data are respectively mapped to d dimensional feature spaces, and the p in the way of Volterra series (p >=1) rank is blocked, so as to obtain N₀(N₀=D-d × τ) individual data frame { (X_t,Y_t) (t=1,2 ..., N₀), and the input as forecast model and target output data, wherein：Input data is X_t=1, x (t), x (t+ τ) ..., x (t+ (d-1) × τ), x²(t),x(t)×x(t+τ),…,x²(t+ (d-1) × τ) ... } target output data be Y_t=x (t+d × τ) (t=1,2 ..., N₀)；

(2), the local optimum positions of primary

Provided with N_SIndividual particle, particle position M dimensional vectors W_iRepresent, population position can use matrixRepresent, i=1,2 ..., N_S；Random initializtion N_SThe initial position W of individual particle_i, and it is Lbest to make each particle local optimum positions_i=W_i；

(3) N, is obtained_SThe initial global optimum position for the colony that individual particle is constituted

I-th (i=1,2 ..., N_S) adaptation value function of the individual particle in searching process be：

${\mathrm{fit}}_i\left(W\right)=\frac{1}{N_0}\underset{t=1}{\overset{N_0}{\Σ}}{\[Y_t-W^T{X}_t\]}^2---\left(1\right)$

The position for choosing the minimum particle of adaptive value is used as N_SThe initial global optimum position Fbest for the colony that individual particle is constituted；

(4), according to universal gravitation algorithm calculate particle i in every one-dimensional space (i=1,2 ..., N_S) accelerationTherefore the acceleration of all particles can use N in M dimension spaces_S× M dimensions matrix a is represented；

(5), the position of more new particle and adaptive value：

$mbest=\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}Fbest\left(i\right)$

$p\left(K\right)=\frac{a_0\·Fbest+a\·fit\left(Fbest\right)}{a_1+a}$

B (K)=1-0.5K/MAXITER

R (K)=[r₁(K),…,r_d(K),…,r_M], (K) wherein $r_d\left(K\right)=\left\{\begin{array}{cc}-1,& if\mathrm{\δ}\left(K\right)>0.5\\ -1,& if\mathrm{\δ}\left(K\right)\≤0.5\end{array}\right.$

$W\left(K\right)=p\left(K\right)+r\left(K\right)b\left(K\right)|mbest-W(K-1)|ln\left(\frac{1}{u\left(K\right)}\right)$

Wherein, δ (K), u (K) are the random number between 0 to 1, a₀,a₁N between respectively 0 to 1_SPeacekeeping N_S× M ties up random vector, and a is the acceleration of particle, and fit (Fbest) is the adaptive value corresponding to global optimum position, i.e. adaptive optimal control value, and K is current iteration step number, and MAXITER is greatest iteration step number；

With particle position W (K) renewal, particle i (i=1,2 ..., N are obtained so as to can be calculated using formula (1)_S) new adaptive value fit (W_i(K))；

(6) local optimum positions of each particle, are updated, if each particle fit (W_i(K)) ＜ fit (Lbest_i), then Lbest_i=W_i(K), otherwise Lbest_iKeep it is constant (i=1,2 ..., N_S)；

Global optimum position is updated, $Fbest=\underset{1\≤i\≤N_S}{\arg min}fit\left({\mathrm{Lbest}}_i\right);$

(7), iterative operation：Judge whether iterative steps reach whether greatest iteration step number MAXITER or adaptive optimal control value fit (Fbest) reach stabilization；If it is, terminating iteration, one group of optimal M dimension core vector is obtainedPerform step (8)；Otherwise, step (4) is returned to continue executing with；

(8) predicted value, is calculated：When predicting the D+1 data, i.e. during prediction data x (D+1), the input data frame for understanding now forecast model according to step (1) is $X_{\hat{t}}={\[x_{\hat{t},1},x_{\hat{t},2},...,x_{\hat{t},M}\]}^T,(\hat{t}=D+1-d\×\mathrm{\τ});$ By the core vector in step (7)Substitute into, so as to obtainMoment, the one-step prediction value of forecast model：

$y_{\hat{t}}=f\left(X_{\hat{t}}\right)=f(x_{\hat{t},1},x_{\hat{t},2},...,x_{\hat{t},M})={\hat{W}}^T{X}_{\hat{t}}.$

Further, it is present invention also offers the method for the acceleration of particle in calculating per the one-dimensional space：

During kth iteration, the inertia mass m of each particle_i(K):

$c_i\left(K\right)=\frac{{\[{\mathrm{fit}}_i\left(W\right)\]}_K-worst\left(K\right)}{best\left(K\right)-worst\left(K\right)}$

$m_i\left(K\right)=\frac{c_i\left(K\right)}{{\mathrm{\Σ}}_{j=1}^{N_S}{c}_i\left(K\right)}$

Wherein, c_i(K) it is the intermediate variable of calculating particle inertia quality, [fit_i(W)]_KParticle i adaptive value when being kth iteration, best (K) and worst (K) represent adaptive optimal control degree and worst fitness；

For the identification of Volterra cores, purpose is vectorial to seek one group of M dimension coreAnd makeMinimum, so in search function minimum value, best (K) and worst (K) are respectively：

According to Formula of Universal Gravitation：F=Gm₁m₂/R², gravitation when obtaining kth iteration between particle i and particle j is：

$F_{ij}^M\left(K\right)=G\left(K\right)\frac{m_i\left(K\right)\×m_j\left(K\right)}{R_{ij}\left(K\right)+\mathrm{\ϵ}}\left(W_j\right(K)-W_i(K\left)\right)$

Wherein, W_i(K),W_j(K) be particle i and particle j in M dimension spaces position, ε is the constant of a very little, R_ij(K) when being kth iteration, particle i and particle j Euclidean distance, i.e. R_ij(K)=| | W_i(K),W_j(K)||₂；G (K) is gravitation coefficient, is represented by：G₀=100, β=20, K are current iteration step number, and MAXITER is greatest iteration step number；

So obtaining the acceleration of particle i in every dimension space

$a_i^M\left(K\right)=\frac{{\mathrm{\Σ}}_{j=1,j\≠i}^{N_S}{h}_j{F}_{ij}^M\left(K\right)}{m_i\left(K\right)}$

Wherein, h_j(1≤j≤N_S, j ≠ i) and represent j-th of random number.

When further, invention further describes predicted value is calculated, multi-step data can also be once predicted, i.e., first prediction timeData, then with the momentPrediction data based on, prediction timeData, and the like, so as to obtain q (q >=1) step advanced prediction value and be：

$y_{\hat{t}+q}=\left\{\begin{array}{cc}f(x_{\hat{t},q+1},...,x_{\hat{t},m},y_{\hat{t}},...,y_{\hat{t}+q-1})& q\∈\{1,2,...,M-1\}\\ f(y_{\hat{t}+q-M},...,y_{\hat{t}+q-1})& q\∈\{M,M+1,...\}\end{array}\right..$

What the goal of the invention of the present invention was realized in：

A kind of Time Series Forecasting Methods based on quantum gravity algorithm of the present invention, build time sequential forecasting models are carried out by using Volterra series expansions.First, on the basis of phase space framing reconfiguration technique is quoted, the input signal vector blocked with Volterra series P ranks is built；Secondly, it is introduced into the quantum gravity algorithm that acceleration is proposed as some variable parameter to MohadesehSoleimanpour et al., the Volterra kernel functions of forecast model is effectively trained using this algorithm；Finally, the input signal vector by the future position moment and the linear combination of Volterra kernel functions, so as to try to achieve the predicted value at the moment.By experimental verification, quantum gravity algorithm is incorporated into the identification of the Volterra kernel functions of nonlinear system, the degree of accuracy of time series forecasting can be effectively improved.

Brief description of the drawings

Fig. 1 is a kind of embodiment flow chart of the Time Series Forecasting Methods of the invention based on quantum gravity algorithm；

Fig. 2 is the position of tri- kinds of sequence particles in searching process of MG, RDD, Laser；

Fig. 3 is learning curve schematic diagram of the present invention under 9 kinds of different pieces of information collection；

Fig. 4 be the present invention to Laser sequences carry out Single-step Prediction predict the outcome and real data comparison diagram；

Fig. 5 be the present invention to Laser sequences carry out 50 step predictions predict the outcome and real data comparison diagram；

When Fig. 6 is different forecast model prediction Laser sequences, the three-dimensional cylinder comparison diagram of its mean square error.

Embodiment

The embodiment to the present invention is described below in conjunction with the accompanying drawings, so that those skilled in the art more fully understands the present invention.Requiring particular attention is that, in the following description, when perhaps the detailed description of known function and design can desalinate the main contents of the present invention, these descriptions will be ignored herein.

Embodiment

Fig. 1 is a kind of embodiment flow chart of the Time Series Forecasting Methods of the invention based on quantum universal gravitation algorithm.

In the present embodiment, as shown in figure 1, the Time Series Forecasting Methods of the invention based on quantum gravity algorithm can be divided into three phases：Data preprocessing phase, forecast model training stage and forecast period.Each stage is described in detail separately below：

(1) data preprocessing phase

S101：If the data set DS={ x (1), x (2) ..., x (D) } of pre-selection, wherein D are data amount check；

First, data set DS smallest embedding dimension number d and delay time T is determined according to the method for differential entropy rate；

Secondly, framing reconstruct is carried out to data set DS using the function windowize in Matlab, and each frame data is respectively mapped to d dimensional feature spaces, and the p in the way of Volterra series (p >=1) rank is blocked, so as to obtain N₀(N₀=D-d × τ) individual data frame { (X_t,Y_t) (t=1,2 ..., N₀), and the input as forecast model and target output data, wherein：Input data is X_t=1, x (t), x (t+ τ) ..., x (t+ (d-1) × τ), x²(t),x(t)×x(t+τ),…,x²(t+ (d-1) × τ) ... } target output data be Y_t=x (t+d × τ) (t=1,2 ..., N₀)；

In the present embodiment, it is the terseness of algorithm, if the dimension size of input data is M, then input data X_t=[x_t,1,x_t,2,…,x_t,M]^T。

(2) the forecast model training stage

S102：Provided with N_SIndividual particle, particle position M dimensional vectors W_iRepresent, population position can use matrixRepresent, i=1,2 ..., N_S；Random initializtion N_SThe initial position W of individual particle_i, and it is Lbest to make each particle local optimum positions_i=W_i。

S103：N is calculated respectively_SThe adaptive value of individual particle, wherein, the i-th (i=1,2 ..., N_S) adaptation value function of the individual particle in searching process be：

${\mathrm{fit}}_i\left(W\right)=\frac{1}{N_0}\underset{i=1}{\overset{N_0}{\Σ}}{\[Y_t-W^T{X}_t\]}^2---\left(1\right)$

The position for choosing the minimum particle of adaptive value is used as N_SThe initial global optimum position Fbest for the colony that individual particle is constituted；

S104：According to universal gravitation algorithm calculate particle i in every one-dimensional space (i=1,2 ..., N_S) accelerationTherefore the acceleration of all particles can use N in M dimension spaces_S× M dimensions matrix a is represented.

Wherein, the acceleration of particle in every one-dimensional space is calculated, detailed process is：

During kth iteration, the inertia mass m of each particle_i(K):

$c_i\left(K\right)=\frac{{\[{\mathrm{fit}}_i\left(W\right)\]}_K-worst\left(K\right)}{best\left(K\right)-worst\left(K\right)}$

$m_i\left(K\right)=\frac{c_i\left(K\right)}{{\mathrm{\Σ}}_{j=1}^{N_S}{c}_i\left(K\right)}$

Wherein, c_i(K) it is the intermediate variable of calculating particle inertia quality, [fit_i(W)]_KParticle i adaptive value when being kth iteration, best (K) and worst (K) represent adaptive optimal control degree and worst fitness；

For the identification of Volterra cores, purpose is vectorial to seek one group of M dimension coreAnd makeMinimum, so in search function minimum value, best (K) and worst (K) are respectively：

$best\left(K\right)=\underset{j\∈\{1,...,N_S\}}{\min }{\[{\mathrm{fit}}_j\left(W\right)\]}_K,worst\left(K\right)=\underset{j\∈\{1,...,N_S\}}{max}{\[{\mathrm{fit}}_j\left(W\right)\]}_K$

According to Formula of Universal Gravitation：F=Gm₁m₂/R², gravitation when obtaining kth iteration between particle i and particle j is：

$F_{ij}^M\left(K\right)=G\left(K\right)\frac{m_i\left(K\right)\×m_j\left(K\right)}{R_{ij}\left(K\right)+\mathrm{\ϵ}}\left(W_j\right(K)-W_i(K\left)\right)$

Wherein, W_i(K),W_j(K) be particle i and particle j in M dimension spaces position, ε is the constant of a very little, R_ij(K) when being kth iteration, particle i and particle j Euclidean distance, i.e. R_ij(K)=| | W_i(K),W_j(K)||₂；G (K) is gravitation coefficient, is represented by：G₀=100, β=20, K are current iteration step number, and MAXITER is greatest iteration step number；

So obtaining the acceleration of particle i in every dimension space

$a_i^M\left(K\right)=\frac{{\mathrm{\Σ}}_{j=1,j\≠i}^{N_S}{h}_j{F}_{ij}^M\left(K\right)}{m_i\left(K\right)}$

Wherein, h_j(1≤j≤N_S, j ≠ i) and represent j-th of random number.

S105:The position of more new particle and adaptive value：

$mbest=\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}Fbest\left(i\right)$

$p\left(K\right)=\frac{a_0\·Fbest+a\·fit\left(Fbest\right)}{a_1+a}$

B (K)=1-0.5K/MAXITER

R (K)=[r₁(K),…,r_d(K),…,r_M], (K) wherein $r_d\left(K\right)=\left\{\begin{array}{cc}-1,& if\mathrm{\δ}\left(K\right)>0.5\\ -1,& if\mathrm{\δ}\left(K\right)\≤0.5\end{array}\right.$

$W\left(K\right)=p\left(K\right)+r\left(K\right)b\left(K\right)|mbest-W(K-1)|ln\left(\frac{1}{u\left(K\right)}\right)$

Wherein, δ (K), u (K) are the random number between 0 to 1, a₀,a₁N between respectively 0 to 1_SPeacekeeping N_S× M ties up random vector, and a is the acceleration of particle, and fit (Fbest) is the adaptive value corresponding to global optimum position, i.e. adaptive optimal control value；K is current iteration step number, and MAXITER is greatest iteration step number；

With particle position W (K) renewal, particle i (i=1,2 ..., N are obtained so as to can be calculated using formula (1)_S) new adaptive value fit (W_i(K))。

S106:The local optimum positions of each particle are updated, if each particle fit (W_i(K)) ＜ fit (Lbest_i), then Lbest_i=W_i(K), otherwise Lbest_iKeep it is constant (i=1,2 ..., N_S)；

Global optimum position is updated, $Fbest=\underset{1\≤i\≤N_S}{\arg min}fit\left({\mathrm{Lbest}}_i\right).$

S107:Iterative operation：Judge whether iterative steps reach whether greatest iteration step number MAXITER or adaptive optimal control value fit (Fbest) reach stabilization；If it is, terminating iteration, one group of optimal M dimension core vector is obtainedPerform step (8)；Otherwise, step (4) is returned to continue executing with.

(3) forecast period

S108：Calculate predicted value：When predicting the D+1 data, i.e. during prediction data x (D+1), the input data frame for understanding now forecast model according to step S101 is $X_{\hat{t}}={\[x_{\hat{t},1},x_{\hat{t},2},...,x_{\hat{t},M}\]}^T,(\hat{t}=D+1-d\×\mathrm{\τ});$ By the core vector in step S107Substitute into, so as to obtainMoment, the one-step prediction value of forecast model：

$y_{\hat{t}}=f\left(X_{\hat{t}}\right)=f(x_{\hat{t},1},x_{\hat{t},2},...,x_{\hat{t},M})={\hat{W}}^T{X}_{\hat{t}}$

In practice, the data of a step once can be only predicted, i.e.,The data at moment, can also once predict multi-step data, i.e. momentData.If once predicting multi-step data, first prediction timeData, then with the momentPrediction data based on, prediction timeData, and the like, so as to obtain q (q >=1) step advanced prediction value and be：

$y_{\hat{t}+q}=\left\{\begin{array}{cc}f(x_{\hat{t},q+1},...,x_{\hat{t},m},y_{\hat{t}},...,y_{\hat{t}+q-1})& q\∈\{1,2,...,M-1\}\\ f(y_{\hat{t}+q-M},...,y_{\hat{t}+q-1})& q\∈\{M,M+1,...\}\end{array}\right..$

Example

In order to illustrate the technique effect of the present invention, the present invention will carry out experimental analysis from emulation and actual chaos time sequence respectively, to weigh the validity of forecast model.Wherein simulation sequence is famous Mackey-Glass (MG) data, several group data sets that actual sequence provides for TimeSeriesPredictionGroup groups:Laser、DailyMinimumTemperatures(DMT)、Electricity_Demand(ED)、CATS_Benchmark(CATS_B)、SunspotNumber(SN)、PolandElectricity(PE)、Dslp.Meanwhile, one group of degeneration measurement data (RDD) that the present invention also utilizes random degradation model to produce is used to verify method proposed by the invention.

For this 9 different emulation and actual chaos time sequence, it is normalized first, then using 2/3 data as training data, 1/3 data are used as test data.In the parameter setting of forecast model, smallest embedding dimension number and time delay are chosen using the method for differential entropy rate；The population number for setting model is 50, and greatest iteration step number MAXITER is 1000, and the exponent number of Volterra series is 3.

Under the conditions of above-mentioned parameter, Fig. 2 (a) (b) (c) sets forth the position of tri- kinds of sequence particles in searching process of MG, RDD, Laser, five-pointed star represents particle optimum position corresponding when adaptive value reaches stable, the optimum position of particle when circle represents each iteration in figure.Although MG sequence randomnesss are weaker, and RDD and Laser sequence randomnesss are stronger, and can be seen that method proposed by the invention from Fig. 2 three width figures is respectively provided with stronger search capability to the random sequence of both types.And then, Fig. 3 gives learning curve schematic diagram of this method under above-mentioned 9 kinds of different pieces of information collection.It can be seen that, in function searching process, for different data sets, this method can possess higher convergence rate and convergence precision.

Training is predicted after terminating using forecast model to test data point.Fig. 4 be using the present invention to Laser sequences carry out Single-step Prediction predict the outcome and real data comparison diagram.Fig. 5 be using the present invention to Laser sequences carry out 50 step predictions predict the outcome and real data comparison diagram.As shown in Figure 4 and Figure 5, real data is represented with solid line, dotted line represents to predict the outcome.As can be seen from Figure 4 and Figure 5, either Single-step Prediction or multi-step prediction, method proposed by the invention all have good estimated performance, and even in the overlength prediction case of 50 steps, it predicts that curve of output also approaches real sequence very much.

Preferably to verify the estimated performance of method proposed by the invention, here single step and multi-step prediction are carried out to above-mentioned 9 different emulation and actual chaos time sequence using mean square error (MSE) and standardization root-mean-square error (NRMSE) respectively as Performance evaluation criterion using forecast model proposed by the invention.In order to ensure the accuracy of test, the prediction to different step-lengths repeats 30 experiments, records the MSE and NRMSE tested every time, its average value is finally sought respectively.

Table 1 is then the error analysis of the model prediction result under different pieces of information collection.

Table 1

The weaker MG sequences of either randomness are can be seen that from the error analysis result of table 1, or for sequences such as randomness very strong DMT, SN, PE, RDD, method proposed by the present invention can obtain preferable prediction effect for the random sequence of both types.

Further, for quantitative measurement and comparison prediction performance, we also use autoregression model (AR), the radial basis neural network (RBF) of standard, optimal beta pruning extreme learning machine model (OPELM), Volterra, cut Lazy learning model (LLpruned), least square method supporting vector machine (LSSVM), and this seven kinds of different forecast models of the invention are predicted to Laser sequences respectively, the training data and test data of all forecast models have all carried out normalized before input model.Intuitively recognize to have to result, when Fig. 6 depicts different forecast model prediction Laser sequences, the three-dimensional cylinder comparison diagram of its mean square error.From comparison diagram as can be seen that model proposed by the invention has higher precision of prediction than popular at present and conventional most of forecast models.

It is above-mentioned test result indicate that, a kind of Time Series Forecasting Methods based on quantum gravity algorithm of the present invention, not only possess higher convergence rate and convergence precision in function searching process, and in the case of single step or multistep advanced prediction, can preferable prediction effect be obtained to different types of random sequence.

Although illustrative embodiment of the invention is described above; in order to which those skilled in the art understand the present invention; it is to be understood that; the invention is not restricted to the scope of embodiment; for those skilled in the art; as long as various change is in the spirit and scope of the present invention that appended claim is limited and is determined, these changes are it will be apparent that all utilize the innovation and creation of present inventive concept in the row of protection.

Claims (3)

1. a kind of Time Series Forecasting Methods based on quantum gravity algorithm, it is characterised in that comprise the following steps：

(1), data prediction

(1.1) the data set DS={ x (1), x (2) ..., x (D) } of pre-selection, is set, wherein D is data amount check；

(1.2) data set DS smallest embedding dimension number d and delay time T, is determined according to the method for differential entropy rate；

(1.3) framing reconstruct, is carried out to data set DS using the function windowize in Matlab, and each frame data are respectively mapped to d dimensional feature spaces, and the p in the way of Volterra series (p >=1) rank is blocked, so as to obtain N₀(N₀=D-d × τ) individual data frame { (X_t,Y_t) (t=1,2 ..., N₀), and the input as forecast model and target output data, wherein：Input data is X_t=1, x (t), x (t+ τ) ..., x (t+ (d-1) × τ), x²(t),x(t)×x(t+τ),…,x²(t+ (d-1) × τ) ... } target output data be Y_t=x (t+d × τ) (t=1,2 ..., N₀)；

(2), the local optimum positions of primary

Provided with N_SIndividual particle, particle position M dimensional vectors W_iRepresent, population position can use matrixRepresent, i=1,2 ..., N_S；The initial position W of the N number of particle of random initializtion_i, and it is Lbest to make each particle local optimum positions_i=W_i；

(3) N, is obtained_SThe initial global optimum position for the colony that individual particle is constituted

I-th (i=1,2 ..., N_S) adaptation value function of the individual particle in searching process be：

${\mathrm{fit}}_i\left(W\right)=\frac{1}{N_0}\underset{t=1}{\overset{N_0}{\Σ}}{\[Y_t-W^T{X}_t\]}^2---\left(1\right)$

The position for choosing the minimum particle of adaptive value is used as N_SThe initial global optimum position Fbest for the colony that individual particle is constituted；

(4), according to universal gravitation algorithm calculate particle i in every one-dimensional space (i=1,2 ..., N_S) accelerationTherefore the acceleration of all particles can use N in M dimension spaces_S× M dimensions matrix a is represented；

(5), the position of more new particle and adaptive value：

$mbest=\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}Fbest\left(i\right)$

$p\left(K\right)=\frac{a_0\·Fbest+a\·fit\left(Fbest\right)}{a_1+a}$

B (K)=1-0.5K/MAXITER

R (K)=[r₁(K),…,r_d(K),…,r_M], (K) wherein $r_d\left(K\right)=\left\{\begin{array}{ccc}-1,& if& \mathrm{\δ}\left(K\right)>0.5\\ -1,& if& \mathrm{\δ}\left(K\right)\≤0.5\end{array}\right.$

$W\left(K\right)=p\left(K\right)+r\left(K\right)b\left(K\right)|mbest-W(K-1)|ln\left(\frac{1}{u\left(K\right)}\right)$

Wherein, δ (K), u (K) are the random number between 0 to 1, a₀,a₁N between respectively 0 to 1_SPeacekeeping N_S× M ties up random vector, and a is the acceleration of particle, and fit (Fbest) is the adaptive value corresponding to global optimum position, i.e. adaptive optimal control value.K is current iteration step number, and MAXITER is greatest iteration step number；

With particle position W (K) renewal, particle i (i=1,2 ..., N are obtained so as to can be calculated using formula (1)_S) new adaptive value fit (W_i(K))；

(6) local optimum positions of each particle, are updated, if each particle fit (W_i(K)) ＜ fit (Lbest_i), then Lbest_i=W_i(K), otherwise Lbest_iKeep it is constant (i=1,2 ..., N_S)；

Global optimum position is updated, $Fbest=\underset{1\≤i\≤N_S}{\arg \min }fit\left({\mathrm{Lbest}}_i\right);$

(7), iterative operation：Judge whether iterative steps reach whether greatest iteration step number MAXITER or adaptive optimal control value fit (Fbest) reach stabilization；If it is, terminating iteration, one group of optimal M dimension core vector is obtainedPerform step (8)；Otherwise, step (4) is returned to continue executing with；

(8) predicted value, is calculated：When predicting the D+1 data, i.e. during prediction data x (D+1), the input data frame for understanding now forecast model according to step (1) is $X_{\hat{t}}={\[x_{\hat{t},1},x_{\hat{t},2},...,x_{\hat{t},M}\]}^T,(\hat{t}=D+1-d\×\mathrm{\τ});$ By the core vector in step (7)Substitute into, so as to obtainMoment, the one-step prediction value of forecast model：

$y_{\hat{t}}=f\left(X_{\hat{t}}\right)=f(x_{\hat{t},1},x_{\hat{t},2},...,x_{\hat{t},M})={\hat{W}}^T{X}_{\hat{t}}.$

2. the Time Series Forecasting Methods according to claim 1 based on quantum gravity algorithm, it is characterised in that in the step (3), calculate the acceleration of particle in every one-dimensional space, detailed process is：

During kth iteration, the inertia mass m of each particle_i(K):

$c_i\left(K\right)=\frac{{\[{\mathrm{fit}}_i\left(W\right)\]}_K-worst\left(K\right)}{best\left(K\right)-worst\left(K\right)}$

$m_i\left(K\right)=\frac{c_i\left(K\right)}{{\Σ}_{j=1}^{N_S}{c}_i\left(K\right)}$

Wherein, c_i(K) it is the intermediate variable of calculating particle inertia quality, [fit_i(W)]_KParticle i adaptive value when being kth iteration, best (K) and worst (K) represent adaptive optimal control degree and worst fitness；

For the identification of Volterra cores, purpose is vectorial to seek one group of M dimension coreAnd makeMinimum, so in search function minimum value, best (K) and worst (K) are respectively：

$best\left(K\right)=\underset{j\∈\left\{1,...,N_S\right\}}{\min }{\[{\mathrm{fit}}_j\left(W\right)\]}_K,worst\left(K\right)=\underset{j\∈\left\{1,...,N_S\right\}}{\max }{\[{\mathrm{fit}}_j\left(W\right)\]}_K$

According to Formula of Universal Gravitation：F=Gm₁m₂/R₂, gravitation when obtaining kth iteration between particle i and particle j is：

$F_{ij}^M\left(K\right)=G\left(K\right)\frac{m_i\left(K\right)\×m_j\left(K\right)}{R_{ij}\left(K\right)+\mathrm{\ϵ}}\left(W_j\right(K)-W_i(K\left)\right)$

Wherein, W_i(K),W_j(K) be grain particle i and particle j in M dimension spaces position, ε is the constant of a very little, R_ij(K) when being kth iteration, particle i and particle j Euclidean distance, i.e. R_ij(K)=| | W_i(K),W_j(K)||₂；G (K) is gravitation coefficient, is represented by：G₀=100, β=20, K are current iteration step number, and MAXITER is greatest iteration step number；

So obtaining the acceleration of particle i in every dimension space

$a_i^M\left(K\right)=\frac{{\Σ}_{j=1,j\≠i}^{N_S}{h}_j{F}_{ij}^M\left(K\right)}{m_i\left(K\right)}$

Wherein, h_j(1≤j≤N_S, j ≠ i) and represent j-th of random number.

3. the Time Series Forecasting Methods according to claim 1 based on quantum gravity algorithm, it is characterised in that in the step (8), multi-step data can also be once predicted when calculating predicted value, i.e., first prediction timeData, then with the momentPrediction data based on, prediction timeData, and the like, so as to obtain q (q >=1) step advanced prediction value and be：

$y_{\hat{t}+q}=\left\{\begin{array}{cc}f\left(x_{\hat{t},q+1},...,x_{\hat{t},m},y_{\hat{t}},...,y_{\hat{t}+q-1}\right)& q\∈\{1,2,...,M-1\}\\ f\left(y_{\hat{t}+q-M},...,y_{\hat{t}+q-1}\right)& q\∈\{M,M+1,...\}\end{array}\right..$