CN104268424A - Comprehensive subway energy consumption forecasting method based on time sequence - Google Patents

Abstract

The invention discloses a comprehensive subway energy consumption forecasting method based on a time sequence. The comprehensive subway energy consumption forecasting method based on the time sequence directly uses sample data of a subway energy management system to build a total subway energy consumption forecasting model based on the time sequence. The comprehensive subway energy consumption forecasting method based on the time sequence gives sufficient consideration to the subway energy consumption periodicity, builds the forecasting model according to the energy consumption forecasting time, and combines a long autoregression model method with a nonlinear least square method to estimate parameters, wherein the long autoregression model method is used for the primary parameter estimation, and the nonlinear least square method is used for the precise parameter estimation. The comprehensive subway energy consumption forecasting method based on the time sequence uses the data recorded in the subway energy management system to directly build the subway energy consumption forecasting model based on the time sequence, does not need to spend energy to research the energy consumption structure of the subway system and avoids the influence to the subway energy consumption forecasting precision due to incomplete consideration to energy consumption influence factors.

Description

A kind of based on seasonal effect in time series metro energy consumption comprehensive prediction method

Technical field

The present invention relates to a kind of Forecasting Methodology of subway energy consumption, belong to energy consumption prediction field, relate to a kind of based on seasonal effect in time series subway energy consumption Forecasting Methodology specifically.

Background technology

Urban track traffic is the important component part of urban public transport system, has freight volume large, and speed is fast, and punctuality rate is high, takes up an area few, pollutes the features such as little, can well solve the congested in traffic problem in current city.Start along with all parts of the country are large batch of metro project, and the operation milimeter number straight line of subway rises, and improve the utilization ratio of the energy for reduction metro operation cost, protection of the environment, energy-saving and emission-reduction are all significant.And in subway power-economizing method, for scientific analysis and the reasonable prediction of energy consumption, not only contribute to optimizing the control strategy of subway and day-to-day operation, also to the energy efficiency management of subway, the energy distribute play a part huge.

Subway energy consumption Forecasting Methodology can be divided into two large classes, and a class is forward model, and a class is data-driven model.Wherein, the former needs there is the model that consumes energy accurately, the energy consumption in train journey of such as subway to each energy consumption equipment of subway; The latter is then from data, sets up effective mathematical model by statistical theory.At present, Chinese scholars and some companies all adopt forward model substantially, have namely done some researchs to the operation energy consumption of subway.Document 1 (Zhang Yanyan. " City Rail Transit System traction and station energy consumption research ") middle author is according to the attribute of train, carry out the modeling and simulating of tractive force energy consumption, and energy consumption calculation has been carried out respectively to energy consumption equipment such as air conditioner, escalator and illuminations in station, establish accurate energy consumption model; Document 2 (Wang Yuming. " quantitative analysis of Rail Transit System energy consumption factor ") author carried out modeling and simulating according to the technical speed of subway, the scheme that stops and load factor, service condition etc. respectively to energy consumption.But, the energy consumption structure of subway system is complicated, and the operation energy consumption of train is a part wherein, in addition also comprises the operation energy consumption at station, here contain a large amount of nonlinear factors, therefore very difficult forward model sets up the comprehensive energy consumption forecast model of subway system.

The energy consumption structure of subway system is complicated, contains a large amount of nonlinear functions, is difficult to the forward model setting up subway total energy consumption.Simultaneously the energy consumption of subway is by the impact of country's legal festivals and holidays, has very large fluctuation, and alternately has obvious mechanical periodicity along with working day and off-day, and the energy consumption under these factors is difficult to express by accurate mathematical model.The present invention is exactly the record data utilizing subway energy management system, proposes a kind of based on seasonal effect in time series subway energy consumption prediction new.

Summary of the invention

For the defective of existing forecast model and the non-linear of subway energy consumption structure itself, the present invention propose a kind of utilize the record data of subway energy management system thus based on seasonal effect in time series statistics set up subway energy consumption forecast model and need not go consider subway system energy consumption structure composition method.

The one that the present invention proposes, based on seasonal effect in time series metro energy consumption comprehensive prediction method, mainly comprises the following steps:

The first step, adopts difference method to subway energy consumption sample sequence x _tadjust, make the sequence after adjusting be stationary time series, be designated as Y;

Second step, adopts Daniel verification, carries out stationarity verification to the stationary time series Y produced in the first step;

3rd step, utilizes second step to be verified the stationary time series Y of the subway energy consumption data obtained by stationarity, and adopt ARMA model to carry out modeling to it, basic representation is:

Wherein { Y _t, t=1,2,3 ... N}, Y _t-1..., Y _t-nfor Y is at t-1 ... the value in t-n moment, ε is average is zero, and variance is σ _ε ²stationary white noise, ε _t, ε _t-1..., ε _t-mfor ε is at t, t-1 ... the value in t-m moment, in formula for autoregressive coefficient, θ _j(j=1,2 ..., m) be running mean coefficient;

4th step, adopt long auto-regression modelling to carry out just estimating to the model parameter of formula 1, then, the essence adopting nonlinear least square method to realize the model parameter of formula 1 is estimated;

Whether the 5th step, utilize the ε in Q criteria tests model to be white noise, if white noise, then illustrate that the model of foundation is applicable to the time series of subway energy consumption sample, if not white noise, then returns the 4th step and re-start parameter estimation;

6th step, utilizes the subway energy consumption model obtained through the 5th step inspection, according to subway energy consumption stationary time series Y t and before the value in moment, to the variable Y in following t+l moment _t+l, (l>0) makes estimation.

Accompanying drawing explanation

Fig. 1 is the building energy consumption prediction process flow diagram of festivals or holidays in the present invention.

Fig. 2 is the BP neural network structure figure calculating energy consumption modified value festivals or holidays in the present invention.

Embodiment

The present invention is based on seasonal effect in time series model to set up and be mainly divided into six steps: comprise model data pre-service, stationary test, model identification and determine the data prediction of rank, the parameter estimation of model, the adaptive test of model and model.Specific implementation is as following:

Step one: the pre-service of model data

Because subway energy consumption data are nonstationary time series, and have certain seasonality and tendency, for raw data, we cannot utilize time series to carry out modeling, need to adjust sample sequence.The present invention adopts Box-Jenkins method, and namely difference method adjusts sample sequence, eliminates its tendency and seasonality, makes the sequence after changing be stationary sequence.

By analyzing subway energy consumption data sequence, we find data with 1 year for period of change, have very strong seasonality.Simultaneously energy consumption is weekly again according to periodically changing on working day (the week), off-day (Saturday, Sunday).Therefore according to the scope of application of model, we have carried out following process to raw data:

If subway energy consumption original data sequence is { x _t(t=1,2 ..., N), the sequence of its to be the cycle be s, can calculus of differences be carried out:

$Y={\▿}_s{X}_t=X_t-X_{t-s},$

Wherein, Y is through differentiated stationary sequence, for s rank difference operator, X _tfor data sequence is in the value of t, X _t-sfor data sequence is in the value in t-s moment.

Time length for subway energy consumption prediction is different, and the present invention sets up following four kinds of forecast models respectively:

(1) hour forecast model (being applicable to the prediction in 24 hours) is set up

$Y={\▿}_{24}{X}_t=X_t-X_{t-24}$

Wherein, Y is through differentiated stationary sequence, be 24 rank difference operators, X _tfor data sequence is in the value of t, X _t-24for data sequence is in the value in t-24 moment.

(2) long-term day forecast model (cycle is 365 days) is set up

$Y={\▿}_{365}{X}_t=(X_t-X_{t-365})$

Wherein, Y is through differentiated stationary sequence, be 365 rank difference operators, X _tfor data sequence is in the value of t, X _t-365for data sequence is in the value in t-365 moment.

(3) short-term day forecast model (cycle is 7 days, have ignored seasonal factor impact) is set up

$Y={\▿}_7{X}_t=X_t-X_{t-7}$

Wherein, Y is through differentiated stationary sequence, be 7 rank difference operators, X _tfor data sequence is in the value of t, X _t-7for data sequence is in the value in t-7 moment.

(4) moon forecast model (being applicable to monthly value prediction) is set up

$Y={\▿}_{12}{X}_t=X_t-X_{t-12}$

Wherein, Y is through differentiated stationary sequence, be 12 rank difference operators, X _tfor data sequence is in the value of t, X _t-12for data sequence is in the value in t-12 moment.

Step 2: stationary test

According to actual needs, select the one in four kinds of forecast models in step one, the stationary sequence Y obtained is carried out to the inspection of stationarity.The method of checking sequence stationarity is a lot, present invention employs Daniel inspection.The Daniel method of inspection is based upon on the basis of Spearman related coefficient.For seasonal effect in time series sample Y, the order of note Y is R _t=R (Y _t), consider that variable is to (t, R _t), t=1,2 ... the Spearman rank correlation coefficient q of N _s, have

wherein N is the number of sample in sequence Y

Structure statistic t obeys t distribution.

Do following test of hypothesis:

The Daniel method of inspection: for level of signifiance α, calculates (t, R by time series Y _t), t=1,2 ... the Spearman rank correlation coefficient q of N _sif, | T|>t _α/2(N-2), then H is refused ₀, think sequence non-stationary, and work as q _sduring >0, think that sequence is on the rise; q _sduring <0, think that sequence has downtrending.When | T|≤t _α/2(N-2), time, H is accepted ₀, can think sequence stationary.Wherein, t _α/2(N-2) value can obtain by looking into t distribution table.

If sequence Y is not steady, then again first difference is carried out to Y, repeat above-mentioned stationary test, until sequence data is stationary sequence.

Step 3: the foundation of model

Utilize the subway day energy consumption stable data sequence Y that step 2 is obtained by stationary test, adopt ARMA model (ARMA, Auto Regressive Moving Average model) to carry out modeling to subway day energy consumption.The basic representation of ARMA (n, m) is:

Wherein { Y _t, t=1,2,3 ... N}, Y _t-1..., Y _t-nfor Y is at t-1 ... the value in t-n moment, residual epsilon is average is zero, and variance is σ _ε ²stationary white noise, ε _t, ε _t-1..., ε _t-mfor ε is at t, t-1 ... the value in t-m moment.In formula for autoregressive coefficient, θ _j(j=1,2 ..., m) be running mean coefficient.

Step 4: the estimation of model parameter

Because the AR model in time series and MA model can regard special arma modeling as, that is: AR is the arma modeling that running mean parameter gets 0, and MA is the arma modeling that auto-regressive parameter gets 0.Can learn according to subway energy consumption actual conditions, the subway energy consumption of a certain day is certainly relevant with proxima luce (prox. luc) or energy consumption data a few days ago, corresponds in time series models, the coefficient of autoregression part certainly be not 0.Therefore, in the present invention, subway energy consumption can not be MA model.No matter adopt AR or arma modeling, the estimation of model parameter can be carried out in the following method.

(1) the first estimation of model parameter

The present invention adopts long auto-regression modelling to realize the first estimation to model parameter, and its process flow diagram is as shown in annex map 1.First estimation comprises following step:

1. long auto-regression model AR (p) is simulated to stationary time series Y,

Z＝φ×W，

Wherein, Z is stationary time series Y value at a time, and W is the value matrix in stationary time series Y other moment before that moment, and φ is its parameter matrix.

Then adopt its parameter matrix of Least Square Method φ, estimated value is denoted as namely for parameter matrix in element, when i is increased to certain numerical value p, all be tending towards 0, then think that the exponent number of autoregression part is p in long auto-regression model, so obtain the parameter matrix of long auto-regression model AR (p) and the inverse function coefficient of long auto-regression model AR (p) $I_i={\widehat{\mathrm{\&phi;}}}_i,(i=\mathrm{1,2},...,p)$ With model order p.

2. get the exponent number that n is ARMA (n, m) model AR part, m is the exponent number of MA part.M=1 is made to start search.

3. because of n+m=p, then n=p-m.Separate following linear system of equations, obtain θ _j(j=1,2 ..., m).

$\left[\begin{array}{c}{I}_{n+1}\\ I_{n+2}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ I_{n+m}\end{array}\right]=\left[\begin{array}{ccccc}{I}_n& I_{n-1}& I_{n-2}& \&CenterDot;\&CenterDot;\&CenterDot;& I_{n+1-m}\\ I_{n+1}& I_n& I_{n-1}& \&CenterDot;\&CenterDot;\&CenterDot;& I_{n+2-m}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ I_{n+m-1}& I_{n+m-2}& I_{n+m-3}& \&CenterDot;\&CenterDot;\&CenterDot;& I_n\end{array}\right]\left[\begin{array}{c}{\mathrm{\&theta;}}_1\\ {\mathrm{\&theta;}}_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&theta;}}_m\end{array}\right]$

4. check | θ _m| whether be less than 10 ^-6, if not, then make m=m+1, return step and 3. circulate; If be less than 10 ^-6, then the θ once circulated before determining _j(j=1,2 ..., m) be running mean parameter, and the order of MA part is m=m-1, then makes n=p-m perform next step.

5. following linear system of equations is separated,

6. check last several value whether be less than 10 ^-6, if so, then save below value, reservation is left for auto-regressive parameter, determine that the order n of AR part equals remaining the number of value, now n<p-m; If be not less than, then need not omit, now n=p-m, for auto-regressive parameter.

Just estimated by above-mentioned, obtain for θ in following essence estimation _jwith initial value

(2) essence of model parameter is estimated

The essence that the present invention adopts nonlinear least square method to carry out implementation model parameter is estimated.

In the iterative computation of nonlinear least square method, first determine the initial value β of iterative computation ⁽⁰⁾and residual epsilon _tinitial value.Wherein, ε _tvalue determined by following formula.

Then carry out essence to model parameter to estimate, its basic thought is, by Taylor series by f (y _t, β) launch, omit the higher order term of more than second order, retain linear term, carry out the Linear least square estimation of iteration, until iteration convergence.The process flow diagram that essence is estimated is with reference to annex map 2, and its concrete calculation procedure is as follows:

1. to stationary time series Y, Confirming model order n is just estimated by parameter, m and model parameter initial value β ⁽⁰⁾and residual epsilon _tinitial value.

2. set k as iterative loop variable, make k=0 start first time iteration.

3. by following various calculating h ^(k), i.e. the difference of model parameter estimation β and its initial value.

For ARMA (n, m) model, have:

Wherein, ε is zero-mean, and variance is σ _ε ²stationary white noise, ε _t-1..., ε _t-mfor ε is at t-1 ... the value in t-m moment, in formula for autoregressive coefficient, θ _j(j=1,2 ..., m) be running mean coefficient.

Make y _t=[Y _t-1y _t-2y _t-nε _t-1ε _t-2ε _t-m] ^t, (m+n=p)

$f(y_t,\mathrm{\&beta;})=y_t^T\mathrm{\&beta;},v_{\mathrm{it}}^{\left(k\right)}=\frac{\&PartialD;f(y_t,\mathrm{\&beta;})}{\&PartialD;{\mathrm{\&beta;}}_i}{|}_{\mathrm{\&beta;}={\mathrm{\&beta;}}^{\left(k\right)}}(i=\mathrm{1,2},...,p)$

$v^{\left(k\right)}=\left[\begin{array}{cccc}{v}_{1,n+1}^{\left(k\right)}& v_{2,n+1}^{\left(k\right)}& \&CenterDot;\&CenterDot;\&CenterDot;& v_{p,n+1}^{\left(k\right)}\\ v_{1,n+2}^{\left(k\right)}& v_{2,n+2}^{\left(k\right)}& \&CenterDot;\&CenterDot;\&CenterDot;& v_{p,n+2}^{\left(k\right)}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ v_{1,N}^{\left(k\right)}& v_{2.N}^{\left(k\right)}& \&CenterDot;\&CenterDot;\&CenterDot;& v_{p,N}^{\left(k\right)}\end{array}\right]$

Then wherein β _ifor i-th component of parameter matrix β, i=1,2 ..., p, β ^(k)for the iterative value in parameter matrix β kth generation, k=0,1 ...

4. h is checked ^(k)whether be less than predetermined precision limits δ (p dimensional vector), if h ^(k)< δ, represents that iteration convergence is in claimed range, then with β=β ^(k)+ h ^(k)as the estimated value of model parameter, iterative computation terminates; If h ^(k)< δ is false, then will carry out next iteration calculating.

5. β is made ^(k+1)=β ^(k)+ h ^(k)as the initial value that next iteration calculates, make k=k+1 proceed to step (3) and continue iteration, until h ^(k)< δ iteration terminates.

Step 5: the adaptive test of model

Complete after model determines rank and parameter estimation, will to test examination to the model set up, whether Basic practice is the ε in testing model is white noise.If white noise, then illustrate that the model of foundation is applicable to sample time-series, if not white noise, then need to re-start and determine rank and parameter estimation.The present invention adopts the Q criterion of model check to carry out model testing.

If the residual error of model of fit is designated as it is ε _testimation.{ ε _tcoefficient of autocorrelation can be calculated by following formula

${\widehat{\mathrm{\&rho;}}}_{\mathrm{\&epsiv;},k}=\frac{\underset{t=k+1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_t{\mathrm{\&epsiv;}}_{t-k}}{\underset{t=k+1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_t^2}(k=\mathrm{1,2},...),$ N is the data length of stationary time series

Mathematical statistics can prove, if { ε _t(t=1,2 ... N) be white noise, when data length N> (200 ~ 300), when k is taken to (20 ~ 30) from 1, be similar to normal distribution, namely approximate have meet standardized normal distribution.Mathematical statistics is pointed out, the quadratic sum of limited standard, separate normal random variable meets χ ²distribution, namely has:

$Q=N\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{\widehat{\mathrm{\&rho;}}}_{\mathrm{\&epsiv;},k}^2,Q~{\mathrm{\&chi;}}^2\left(l\right)$

But due to { ε _tcalculated by model, receive the constraint of n+m model parameter, thus χ ²the degree of freedom of variable is l-m-n, generally for the sake of simplicity, gets l-m-n=30.When getting level of confidence and being 95%, by χ ²table checks in χ ²(30) ≈ 44, forms inspection formula Q<=44 like this, when Q value meets this formula, then thinks that corresponding model is applicable models.

Step 6: the data prediction of model

Established the arma modeling of applicable sample time-series by above-mentioned steps, next will forecast time series.It is according to stationary time series Y t and before the value i.e. (Y in moment _t, Y _t-1...) variable Y to the following t+l moment _t+l, (l>0) makes estimation, and estimator is denoted as it is Y _t, Y _t-1... linear combination.As follows time series is forecast.

For ARMA (n, m) model, reverse form according to it obtain the forecast of l step

Formula is: ${\hat{Y}}_t\left(l\right)=\underset{j=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{F_j}^{\left(l\right)}{Y}_{t+1-j}$

Wherein, Y _t+1-jfor stationary time series Y is in t+1-j moment value, coefficient F _j ^(l)can by inverse function { F _j, j=1,2 ... determine, computing formula is as follows:

Inverse function computing formula: θ and be respectively the coefficient in arma modeling

Coefficient F _j ^(l)computing formula: ${F_j}^{\left(l\right)}=\left\{\begin{array}{c}{F}_j,l=1\\ F_{j+l-1}+\underset{i=1}{\overset{l-1}{\mathrm{\&Sigma;}}}{F}_i{F_j}^{(l-i)},l\&NotEqual;1\end{array}\right.$

Beneficial effect of the present invention: by utilizing the record data of subway energy management system, based on seasonal effect in time series analytical approach, directly set up subway energy consumption forecast model, do not need to require efforts and study the energy consumption structure of subway system, more can not affect the energy consumption precision of prediction of subway because of the inconsiderate complete of energy consumption factor.

Claims (4)

1., based on a seasonal effect in time series metro energy consumption comprehensive prediction method, mainly comprise the following steps:

The first step, adopts difference method to subway energy consumption sample sequence x _tadjust, make the sequence after adjusting be stationary time series, be designated as Y;

Second step, adopts Daniel verification, carries out stationarity verification to the stationary time series Y produced in the first step;

3rd step, utilizes second step to be verified the stationary time series Y of the subway energy consumption data obtained by stationarity, and adopt ARMA model to carry out modeling to it, basic representation is:

(formula 1)

Wherein { Y _t, t=1,2,3 ... N}, Y _t-1..., Y _t-nfor Y is at t-1 ... the value in t-n moment, ε is average is zero, and variance is σ _ε ²stationary white noise, ε _t, ε _t-1..., ε _t-mfor ε is at t, t-1 ... the value in t-m moment, in formula for autoregressive coefficient, θ _j(j=1,2 ..., m) be running mean coefficient;

4th step, adopt long auto-regression modelling to carry out just estimating to the model parameter of formula 1, then, the essence adopting nonlinear least square method to realize the model parameter of formula 1 is estimated;

Whether the 5th step, utilize the ε in Q criteria tests model to be white noise, if white noise, then illustrate that the model of foundation is applicable to the time series of subway energy consumption sample, if not white noise, then returns the 4th step and re-start parameter estimation;

6th step, utilizes the subway energy consumption model obtained through the 5th step inspection, according to subway energy consumption stationary time series Y t and before the value in moment, to the variable Y in following t+l moment _t+l, (l>0) makes estimation.

2. metro energy consumption comprehensive prediction method as claimed in claim 1, wherein in a first step, if subway energy consumption sample original series is { x _t(t=1,2 ..., N), and its cycle be s, can calculus of differences be carried out:

$Y={\&dtri;}_s{X}_t=X_t-X_{t-s},$

Wherein, Y is through differentiated stationary sequence, for s rank difference operator, X _tfor data sequence is in the value of t, X _t-sfor data sequence is in the value in t-s moment.

3. metro energy consumption comprehensive prediction method as claimed in claim 1, when wherein carrying out just estimating in the 4th step,

Step one, simulates long auto-regression model AR (p) to stationary time series Y,

Z＝φ×W，

Wherein, Z is stationary time series Y value at a time, and W is the value matrix in stationary time series Y other moment before that moment, and φ is its parameter matrix;

Then adopt its parameter matrix of Least Square Method φ, estimated value is denoted as namely for parameter matrix in element, when i is increased to certain numerical value p, all be tending towards 0, then think that the exponent number of autoregression part is p in long auto-regression model, so the model parameter of obtaining (i.e. inverse function) $I_i=\widehat{{\mathrm{\&phi;}}_i},(i=\mathrm{1,2},...,p)$ With model order p;

Step 2, gets the exponent number that n is ARMA (n, m) model AR part, and m is the exponent number of MA part; M=1 is made to start search;

Step 3, because of n+m=p, then n=p – m, separates following linear system of equations, obtains θ _j(j=1,2 ..., m):

$\left[\begin{array}{c}{I}_{n+1}\\ I_{n+2}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ I_{n+m}\end{array}\right]=\left[\begin{array}{ccccc}{I}_n& I_{n-1}& I_{n-2}& \&CenterDot;\&CenterDot;\&CenterDot;& I_{n+1-m}\\ I_{n+1}& I_n& I_{n-1}& \&CenterDot;\&CenterDot;\&CenterDot;& I_{n+2-m}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ I_{n+m-1}& I_{n+m-2}& I_{n+m-3}& \&CenterDot;\&CenterDot;\&CenterDot;& I_n\end{array}\right]\left[\begin{array}{c}{\mathrm{\&theta;}}_1\\ {\mathrm{\&theta;}}_2\\ \&CenterDot;\&CenterDot;\&CenterDot;\\ {\mathrm{\&theta;}}_m\end{array}\right]$

Step 4, checks | θ _m| whether level off to 0, if not, then make m=m+1, return step 3 circulation; If be tending towards 0, then the θ once circulated before determining _j(j=1,2 ..., m) be running mean parameter, and the order of MA part is m=m-1, then makes n=p-m perform next step;

Step 5, separates following linear system of equations,

Step 6, checks last several value whether be less than 10 ^-6, if so, then save below value, reservation is left for auto-regressive parameter, determine that the order n of AR part equals remaining the number of value, now n<p-m; If be not less than, then need not omit, now n=p-m, (i=1,2 ..., n) be auto-regressive parameter.

4. metro energy consumption comprehensive prediction method as claimed in claim 1, wherein in the 6th step,

For ARMA (n, m) model, reverse form according to it obtaining l step prediction formula is:

${\hat{Y}}_i\left(l\right)=\underset{j=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{F_j}^{\left(l\right)}{Y}_{t+1-j}$

Wherein, Y _t+1-jfor stationary sequence Y is in t+1-j moment value, coefficient can by inverse function { F _j, j=1,2 ... determine, computing formula is as follows:

Inverse function computing formula: θ and be respectively the coefficient in arma modeling

Coefficient F _j ^(l)computing formula: ${F_j}^{\left(l\right)}=\left\{\begin{array}{c}{F}_j,l=1\\ F_{j+l-1}+\underset{i=1}{\overset{l-1}{\mathrm{\&Sigma;}}}{F}_i{F_j}^{(l-i)},l\&NotEqual;1\end{array}\right..$