CN104699991A - Grey system theory-based method of predicting annual heat supply of urban heating system - Google Patents

Abstract

The invention relates to a method of predicting annual heat supply of an urban heating system, in particular to a grey system theory-based method of predicting annual heat supply of the urban heating system, belongs to the technical field of load prediction for heating systems and solves the problem that the existing prediction method of the single heating system is unable to be used in the planning of urban heating systems. The method is technically characterized by including: selecting a historic value of an initial variable and a historic value of a dependent variable; checking relational degree between the initial variable and the dependent variable; subjecting the initial variable to collinearity test, and performing screening to obtain an explanatory variable; building a GM (1, 1) and GM (1, N+1) model based on the grey system theory, and solving parameters of the model; building a grey system state equation to predict annual heat supply in a planned period. The method is applicable to the planning of urban energy systems or heating systems.

Description

Based on the urban heating system year heating load Forecasting Methodology of gray system theory

Technical field

The present invention relates to a kind of urban heating system year heating load Forecasting Methodology, particularly relate to a kind of urban heating system year heating load Forecasting Methodology based on gray system theory, belong to heating system load prediction technical field.

Background technology

Along with the raising of China's Development of China's Urbanization and living standards of the people, the comfort level of Building Indoor Environment and controllability are had higher requirement, China's urban heating area and heating system heating load constantly increase, nineteen ninety to 2012 year, China's Areas benefiting from central heating rise to 51.84 hundred million square metres from 2.13 hundred million square metres, year heat supply total amount is increased to 29.5 hundred million lucky Jiao from 2.87 hundred million lucky Jiao, and heating system energy consumption has become the important component part of social total energy consumption.But China's anxiety energy supply present situation and greatly develop low-carbon economy background under, energy structure of need making rational planning for, improve energy efficiency.

On the basis that urban energy system planning is based upon energy demand reasonable prediction, heat supply planning is mainly divided into two classes: a class is for single heating system, and based on enlarging scale recent and at a specified future date, in prediction, long-term heating demand, planning system extends capacity; Another kind of is for cities and towns, in original dispersion or central heating basis, consider cities and towns in the recent period and long term planning development, prediction heating demand or year heating load, the Energy output of whole heating system and supply in planning cities and towns.The domestic and international load prediction for single heating system is at present often based on the heating power of objective system and the historical data such as operational factor, out door climatic parameter of hydraulic regime, the prediction of urban heating system loading cannot be used for, and when urban heating systems organization, main employing index budgetary estimate method or Degree Days prediction thermal load or year heating load, estimation basis source is single, can not reflect the development leveies such as the economy in cities and towns, population, industry and construction.

Therefore, need one badly and can be applicable to cities and towns energy resource system or heating system planning, and with town planning according to consistent, Heat Supply Engineering development level can not only be embodied, also can reflect the heating system year heating load Forecasting Methodology of the development leveies such as cities and towns economy, population, industry and construction.

Summary of the invention

The object of the invention is to propose a kind of urban heating system year heating load Forecasting Methodology based on gray system theory, the problem of urban heating systems organization can not be used for the Forecasting Methodology solved for single heating system.

The present invention for solving the problems of the technologies described above adopted technical scheme is:

Urban heating system year heating load Forecasting Methodology based on gray system theory of the present invention, realizes according to following steps:

Step one, selection initializaing variable, the numerical value over the years of dependent variable: statistics cities and towns year over the years heating load, as dependent variable, adds up town planning index of correlation over the years as initializaing variable;

The degree of association of step 2, the selected initializaing variable of inspection and dependent variable: apply the grey relational grade that grey Relational Analysis Method calculates each initializaing variable and dependent variable, checked the degree of association of selected initializaing variable and dependent variable by the grey relational grade calculated, preliminary screening is carried out to initializaing variable simultaneously;

Step 3, initializaing variable collinearity are checked, and filter out explanatory variable: adopt correlation coefficient matrix inspection initializaing variable collinearity, then by stepwise regression analysis, screening draws explanatory variable;

Step 4, build GM (1,1) based on gray system theory and GM (1, N+1) model, solving model parameter;

Step 5, set up gray system state equation, the year heating load of prediction project period.

The invention has the beneficial effects as follows:

1. the explanatory variable chosen of the present invention and heating system operational factor have nothing to do, not by the restriction of heating system operational factor in cities and towns, and closely related with town development situation, the development leveies such as urban heating engineering, economy, population, industry and construction can be reflected, with town planning according to consistent.

2. the present invention take cities and towns as heating range, not by the restricted number of heating system in cities and towns, the year heating load of prediction reflects whole urban heating load, simultaneously according to meteorologic parameter, can be the thermal load under heating outdoor calculate temperature by year heating load conversion, as the foundation of cities and towns thermal source planning.

3. the present invention eliminates correlativity between initializaing variable by screening, simplifies as some explanatory variables, can Simplified prediction model, and what ensure dependent variable is explanatory simultaneously.

4. the present invention is based on gray system theory set up state equation be suitable for sample data less cities and towns year heating load prediction, the present invention simultaneously possesses study, can according to Data Update year by year, again explanatory variable is screened and computation model parameter, set up new state equation, ensure real-time and the reliability of prediction.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of the inventive method;

Fig. 2 is Harbin City 2003-2012 heating load statistical value and 2011-2020 heating load predicted value schematic diagram in the embodiment of the present invention.

Embodiment

Further describe the specific embodiment of the present invention by reference to the accompanying drawings.

Embodiment one: present embodiment is described below in conjunction with Fig. 1, a kind of urban heating system year heating load Forecasting Methodology based on gray system theory described in present embodiment, is characterized in that said method comprising the steps of:

Step one, selection initializaing variable, the numerical value over the years of dependent variable: with urban heating system year heating load for forecasting object, statistics cities and towns year over the years heating load is as dependent variable, because heating system is as a part for northern cities and towns infrastructure construction, except the development level of Heat Supply Engineering self, also be subject to the restriction of the development leveies such as cities and towns economy, population, industry and construction, therefore add up town planning index of correlation over the years as initializaing variable;

Described index of correlation is: heat capacity, area of heat-supply service, heat supply pipeline length, residential estate, public facilities, industrial land, city total population in the end of the year, the total amount in city, regional GDP, city be regional GDP, city primary industry total output value, city secondary industry total output value, city tertiary industry total output value, residential households per capita disposable income, town dweller's per capita consumption expenditure per capita, urban house building operation area and urban house building completed floor space, the statistical value of totally 17 initializaing variables.

The degree of association of step 2, the selected initializaing variable of inspection and dependent variable: apply the grey relational grade that grey Relational Analysis Method calculates each initializaing variable and dependent variable, the degree of association of selected initializaing variable and dependent variable is checked by the grey relational grade calculated, preliminary screening is carried out to initializaing variable simultaneously, reject the initializaing variable that grey relational grade is lower;

Step 3, initializaing variable collinearity are checked, and filter out explanatory variable: because the collinearity that may exist between initializaing variable understands the accuracy of impact prediction model.Therefore adopt correlation coefficient matrix inspection initializaing variable collinearity, then by stepwise regression analysis, under the prerequisite taking into account and model accuracy explanatory to dependent variable, screening draws explanatory variable;

Step 4, build GM (1,1) based on gray system theory and GM (1, N+1) model, solving model parameter;

Step 5, set up gray system state equation, the year heating load of prediction project period.

Embodiment two: present embodiment and embodiment one unlike: the computation process of the grey relational grade described in step 2 is as follows:

Step 2 one, employing averaging method carry out nondimensionalization to dependent variable and each initializaing variable:

$X_{i,0}=\frac{x_{i,0}}{{\stackrel{\‾}{x}}_0},X_{i,j}=\frac{x_{i,j}}{{\stackrel{\‾}{x}}_j},(i=\mathrm{1,2},...,m;j=\mathrm{1,2},...,n)---\left(1\right)$

Wherein, X _{i, 0}, x _{i, 0}be respectively dimensionless number and the original value of the i-th annual dependent variable, X _{i, j}, x _{i, j}be respectively dimensionless number and the original value of an i-th annual jth initializaing variable, n is initializaing variable number, and m is the statistics year number of degrees, the average in dependent variable and a jth initializaing variable m year respectively;

Step 2 two, with dimensionless dependent variable for auxiliary sequence X ₀=(X _1,0, X _2,0..., X _{m, 0}) ^t, with dimensionless initializaing variable for subsequence X _j=(X _{1, j}, X _{2, j}..., X _{m, j}) ^t, calculate each annual auxiliary sequence X ₀with each subsequence X _jabsolute value difference δ _{i, j}, and get maximal value δ wherein _maxwith minimum value δ _min:

δ _i，j＝|X _i，0-X _i，j| (2)

δ _max＝max{δ _i，j；i＝1，2，…，m；j＝1，2，…，n} (3)

δ _min＝min{δ _i，j；i＝1，2，…，m；j＝1，2，…，n} (4)

Step 2 three, calculate each annual auxiliary sequence X ₀with each subsequence X _jgrey incidence coefficient l _{i, j}:

$l_{i,j}=\frac{{\mathrm{\δ}}_{\max }+{\mathrm{\ζ}\×\mathrm{\δ}}_{\min }}{{\mathrm{\δ}}_{i,j}+{\mathrm{\ζ}\×\mathrm{\δ}}_{\max }}---\left(5\right)$

Wherein, ζ is resolution ratio, and its value is 0.5;

The grey relational grade of step 2 four, each initializaing variable and dependent variable is:

$L_j=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{l}_{i,j}---\left(6\right).$

Embodiment three: present embodiment and embodiment one or two unlike:

Described in step 2 to carry out the process of preliminary screening to initializaing variable as follows:

As grey relational grade L _jbe greater than grey relational grade and judge limit value L ^*time, judge that the degree of association of a jth initializaing variable and dependent variable is high, otherwise reject this initializaing variable, complete preliminary screening, wherein L ^*value is 0.6 ~ 0.8.

Embodiment four: one of present embodiment and embodiment one to three unlike: the detailed process of initializaing variable collinearity described in step 3 inspection is as follows:

Calculate the relative coefficient between initializaing variable:

$r_{\mathrm{uv}}=\frac{\mathrm{Cov}(X_u,X_v)}{\sqrt{\mathrm{Var}\left(X_u\right)}\·\sqrt{\mathrm{Var}\left(X_v\right)}},(u=\mathrm{1,2},...,n;v=\mathrm{1,2},...,n)---\left(7\right)$

Wherein, r _uvbe the relative coefficient between u and v initializaing variable, Var (X _u), Var (X _v) be respectively subsequence X _u, X _vvariance, Cov (X _u, X _v) be subsequence X _u, X _vcovariance;

Relative coefficient between all initializaing variables forms symmetrical correlation coefficient matrix r, and diagonal entry value is 1, more than matrix r diagonal line or the element absolute value of certain row or column lower be designated as | r _pq|, when it is greater than correlation prediction limit value r ^*time, judge u initializaing variable and v initializaing variable correlativity strong, when having strong correlation between at least one initializaing variable and other initializaing variables, judge to there is multicollinearity between initializaing variable, then adopt the collinearity between Stepwise Regression Method elimination initializaing variable, filter out explanatory variable, otherwise all initializaing variables are explanatory variable, wherein r simultaneously ^*value 0.8.

Embodiment five: one of present embodiment and embodiment one to four unlike: the detailed process filtering out explanatory variable described in step 3 is as follows:

Stepwise Regression Method is adopted to screen initializaing variable: the horizontal α of specification test, if the probable value p that in regression model, the t test statistics of estimates of parameters is corresponding is less than insolation level α, then retain this initializaing variable, otherwise delete this initializaing variable, the t inspection of all initializaing variables is completed according to forward direction or backward stepwise regression analysis order, the initializaing variable finally retained is the explanatory variable of screening and obtaining, and is designated as z ₁, z ₂..., z _n, wherein α value is 5% ~ 10%.

Embodiment six: one of present embodiment and embodiment one to five unlike: the structure described in step 4 is based on GM (1, the 1) model of gray system theory, and the detailed process of solving model parameter is:

First pre-service is carried out to statistics, the original series of the calendar year statistics value of dependent variable and each explanatory variable composition is designated as with wherein k=1,2 ..., N, carries out the one-accumulate formation sequence obtaining dependent variable and each explanatory variable after adding up, is designated as with

$x_0^{\left(1\right)}=(x_{\mathrm{1,0}},\underset{s=1}{\overset{2}{\mathrm{\Σ}}}{x}_{s,0},...,\underset{s=1}{\overset{i}{\mathrm{\Σ}}}{x}_{s,0},...\underset{s=1}{\overset{m}{\mathrm{\Σ}}}{x}_{s,0})---\left(8\right)$

$z_k^{\left(1\right)}=(z_{1,k},\underset{s=1}{\overset{2}{\mathrm{\Σ}}}{z}_{s,k},...,\underset{s=1}{\overset{i}{\mathrm{\Σ}}}{z}_{s,k},...,\underset{s=1}{\overset{m}{\mathrm{\Σ}}}{z}_{s,k})---\left(9\right)$

Then right set up GM (1,1) model respectively, obtain GM (1, the 1) model of each explanatory variable:

$\frac{{\mathrm{dz}}_k^{\left(1\right)}}{\mathrm{d\τ}}=-a_{\mathrm{kk}}{z}_k^{\left(1\right)}+c_k---\left(10\right)$

Wherein, a _kk, c _kfor the parameter of a kth explanatory variable GM (1,1) model, k=1,2 ..., N, τ are the time;

Discrete differential equation, obtaining matrix equation is:

Y _k＝B _kA _k(11)

Wherein, $Y_k=\left[\begin{array}{c}{z}_{2,k}^{\left(0\right)}\\ z_{3,k}^{\left(0\right)}\\ .\\ .\\ .\\ z_{m,k}^{\left(0\right)}\end{array}\right],B_k=\left[\begin{array}{cc}-\frac{1}{2}[z_{1,k}^{\left(1\right)}+z_{2,k}^{\left(1\right)}]& 1\\ -\frac{1}{2}[z_{2,k}^{\left(1\right)}+z_{3,k}^{\left(1\right)}]& 1\\ .& .\\ .& .\\ .& .\\ -\frac{1}{2}[z_{(m-1),k}^{\left(1\right)}+z_{m,k}^{\left(1\right)}]& 1\end{array}\right],A_k=\left[\begin{array}{c}{a}_{\mathrm{kk}}\\ c_k\end{array}\right];$

GM (1, the 1) model parameter of a kth explanatory variable is solved according to least square method:

${\hat{A}}_k=\left[\begin{array}{c}{\hat{a}}_{\mathrm{kk}}\\ {\hat{c}}_k\end{array}\right]={\left(B_k^T{B}_k\right)}^{-1}{B}_k^T{Y}_k---\left(12\right)$

Wherein, for A _k, a _kk, c _kestimated value, k=1,2 ..., N.

Embodiment seven: one of present embodiment and embodiment one to six unlike: the structure described in step 4 is based on GM (1, the N+1) model of gray system theory, and the detailed process of solving model parameter is:

To the one-accumulate formation sequence of dependent variable set up GM (1, the N+1) model about N number of explanatory variable one-accumulate formation sequence, obtain GM (1, the N+1) model of dependent variable:

GM (1, N+1) model reflects the impact for dependent variable first order derivative of N number of explanatory variable and dependent variable itself.Its differential equation form is:

$\frac{{\mathrm{dx}}_0^{\left(1\right)}}{\mathrm{d\τ}}=-a_{00}{x}_0^{\left(1\right)}+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{a}_{0k}{z}_k^{\left(1\right)}---\left(13\right)$

Wherein, a ₀₀, a _0kfor the parameter of GM (1, N+1) model, k=1,2 ..., N;

Discrete differential equation, obtaining matrix equation is:

Y ₀＝B ₀A ₀(14)

Wherein,

$Y_0=\left[\begin{array}{c}{x}_{\mathrm{2,0}}^{\left(0\right)}\\ x_{\mathrm{3,0}}^{\left(0\right)}\\ .\\ .\\ .\\ x_{m,0}^{\left(0\right)}\end{array}\right..,\left[\begin{array}{cccc}-\frac{1}{2}[x_{\mathrm{1,0}}^{\left(1\right)}+x_{\mathrm{2,0}}^{\left(1\right)}]& z_{\mathrm{2,1}}^{\left(1\right)}& ...& z_{2,N}^{\left(1\right)}\\ -\frac{1}{2}[x_{\mathrm{2,0}}^{\left(1\right)}+x_{\mathrm{3,0}}^{\left(1\right)}]& z_{\mathrm{3,1}}^{\left(1\right)}& ...& z_{3,N}^{\left(1\right)}\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ -\frac{1}{2}[x_{(m-1),0}^{\left(1\right)}+x_{m,0}^{\left(1\right)}]& z_{m,1}^{\left(1\right)}& ...& z_{m,N}^{\left(1\right)}\end{array}\right],A_0=\left[\begin{array}{c}{a}_{00}\\ a_{01}\\ .\\ .\\ .\\ a_{0N}\end{array}\right];$

GM (1, N+1) model parameter is solved according to least square method:

${\hat{A}}_0=\left[\begin{array}{c}{\hat{a}}_{00}\\ {\hat{a}}_{01}\\ .\\ .\\ .\\ {\hat{a}}_{0N}\end{array}\right]={\left(B_0^T{B}_0\right)}^{-1}{B}_0^T{Y}_0---\left(15\right)$

Wherein, for A ₀, a ₀₀, a _0kestimated value, k=1,2 ..., N.

Embodiment eight: one of present embodiment and embodiment one to seven unlike: set up gray system state equation described in step 5, prediction project period year heating load detailed process be:

The GM (1 of the N number of explanatory variable of simultaneous, 1) model (formula (10)), the GM (1 of dependent variable, N+1) model (formula (13)), formula (12) and formula (15) are substituted into respectively in formula (10) and formula (13), the state equation obtaining gray system is:

The calendar year statistics value of dependent variable and explanatory variable is initial value, the setting year project period number of degrees, and application Runge-Kutta method solving state equation obtains the one-accumulate formation sequence of dependent variable carry out the predicted value that regressive reduction obtains dependent variable again

${\hat{x}}_0=({\hat{x}}_{\mathrm{1,0}}^{\left(1\right)},{\hat{x}}_{\mathrm{2,0}}^{\left(1\right)}-{\hat{x}}_{\mathrm{1,0}}^{\left(1\right)},...,{\hat{x}}_{i,0}^{\left(1\right)}-\underset{s=1}{\overset{i-1}{\mathrm{\Σ}}}{\hat{x}}_{s,0}^{\left(1\right)},...,{\hat{x}}_{M,0}^{\left(1\right)}-\underset{s=1}{\overset{M-1}{\mathrm{\Σ}}}{\hat{x}}_{s,0}^{\left(1\right)})---\left(17\right)$

Wherein, M is the total year number of degrees of statistics and project period.

The embodiment of the present invention is as follows:

Illustrate further below in conjunction with the prediction of Harbin City's year heating load.

For verifying the precision of Forecasting Methodology of the present invention, get the year heating load of 2003 ~ 2010 years statistics prediction 2011 ~ the year two thousand twenties, wherein heating load predicted values in 2011 and 2012 and statistical value contrast, computational prediction error.

Step one, selection initializaing variable, the numerical value over the years of dependent variable.

Harbin City 2003-2012 heating load (x is taken passages according to " Chinese town and country construction database " ₀) as dependent variable statistical value, in table 1.

Table 1 Harbin City 2003-2012 heating load statistical value

Harbin City 2003-2010 heat capacity (x is taken passages according to " Urban Construction in China statistical yearbook ", " regional economic statistics yearbook " and " Chinese town and country construction database " ₁), area of heat-supply service (x ₂), heat supply pipeline length (x ₃), residential estate (x ₄), public facilities (x ₅), industrial land (x ₆), city total population in the end of the year (x ₇), the total amount (x in city ₈), regional GDP (x ₉), city regional GDP (x per capita ₁₀), city primary industry total output value (x ₁₁), city secondary industry total output value (x ₁₂), city tertiary industry total output value (x ₁₃), residential households per capita disposable income (x ₁₄), town dweller's per capita consumption expenditure (x ₁₅), urban house building operation area (x ₁₆) and urban house building completed floor space (x ₁₇) etc. the statistical value of 17 initializaing variables, in table 2.

Table 2 Harbin City 2003-2010 initializaing variable statistical value

The degree of association of step 2, the selected initializaing variable of inspection and dependent variable.

Adopt averaging method to dependent variable x ₀with initializaing variable x ₁~ x ₁₇carry out nondimensionalization, form auxiliary sequence X ₀with subsequence X _j, calculate the absolute value difference of 2003 ~ 2010 annual each subsequences and auxiliary sequence, obtain its maximal value δ _max=0.40862, minimum value δ _min=0.00146.

Resolution ratio ζ gets 0.5, calculates 2003 ~ 2010 annual auxiliary sequence X ₀with each subsequence X _jgrey incidence coefficient l _{i, j}, obtain the grey relational grade L of each initializaing variable and dependent variable _j, in table 3.

Table 3 initializaing variable x ₁~ x ₁₇with dependent variable x ₀grey relational grade L ₁~ L ₁₇

Grey relational grade	L 1	L 2	L 3	L 4	L 5	L 6	L 7	L 8	L 9
										Numerical value	0.9437	0.9519	0.8997	0.8357	0.8087	0.8072	0.7274	0.7657	0.8301
Grey relational grade	L 10	L 11	L 12	L 13	L 14	L 15	L 16	L 17
										Numerical value	0.8456	0.9170	0.8280	0.7896	0.9013	0.8931	0.6089	0.6762

Grey relational grade judges limit value L ^*get 0.6, the degree of association L of each initializaing variable and dependent variable _jall higher than 0.6, illustrate selected 17 initializaing variables and the dependent variable degree of association high.

Step 3, initializaing variable collinearity are checked, and filter out explanatory variable.

Calculate relative coefficient between initializaing variable, form correlation coefficient matrix r, in table 4.

The correlation coefficient matrix r of table 4 initializaing variable

	x 1	x 2	x 3	x 4	x 5	x 6	x 7	x 8	x 9	x 10	x 11	x 12	x 13	x 14	x 15	x 16	x 17
																		x 1	1.00	1.00	1.00	0.91	0.86	0.92	0.97	0.98	0.98	0.98	0.98	0.97	0.98	0.99	0.97	0.53	0.30
x 2	1.00	1.00	1.00	0.92	0.87	0.93	0.97	0.99	0.99	0.99	0.99	0.98	0.99	0.99	0.98	0.54	0.33

x 3	1.00	1.00	1.00	0.90	0.84	0.91	0.96	0.99	0.98	0.98	0.97	0.98	0.98	0.99	0.97	0.53	0.31
																		x 4	0.91	0.92	0.90	1.00	0.99	1.00	0.97	0.87	0.86	0.86	0.94	0.86	0.85	0.87	0.84	0.30	0.09
x 5	0.86	0.87	0.84	0.99	1.00	0.99	0.94	0.81	0.80	0.80	0.90	0.80	0.79	0.81	0.78	0.24	0.04
																		x 6	0.92	0.93	0.91	1.00	0.99	1.00	0.97	0.89	0.87	0.87	0.94	0.86	0.86	0.88	0.85	0.32	0.09
x 7	0.97	0.97	0.96	0.97	0.94	0.97	1.00	0.93	0.94	0.93	0.99	0.93	0.93	0.94	0.91	0.43	0.20
																		x 8	0.98	0.99	0.99	0.87	0.81	0.89	0.93	1.00	0.99	0.99	0.96	0.98	0.99	0.99	0.99	0.59	0.40
x 9	0.98	0.99	0.98	0.86	0.80	0.87	0.94	0.99	1.00	1.00	0.97	1.00	1.00	1.00	0.99	0.62	0.42
																		x 10	0.98	0.99	0.98	0.86	0.80	0.87	0.93	0.99	1.00	1.00	0.97	1.00	1.00	1.00	0.99	0.62	0.43
x 11	0.98	0.99	0.97	0.94	0.90	0.94	0.99	0.96	0.97	0.97	1.00	0.97	0.97	0.98	0.96	0.55	0.34
																		x 12	0.97	0.98	0.98	0.86	0.80	0.86	0.93	0.98	1.00	1.00	0.97	1.00	0.99	0.99	0.99	0.60	0.39
x 13	0.98	0.99	0.98	0.85	0.79	0.86	0.93	0.99	1.00	1.00	0.97	0.99	1.00	1.00	1.00	0.64	0.46
																		x 14	0.99	0.99	0.99	0.87	0.81	0.88	0.94	0.99	1.00	1.00	0.98	0.99	1.00	1.00	0.99	0.62	0.42
x 15	0.97	0.98	0.97	0.84	0.78	0.85	0.91	0.99	0.99	0.99	0.96	0.99	1.00	0.99	1.00	0.66	0.50
																		x 16	0.53	0.54	0.53	0.30	0.24	0.32	0.43	0.59	0.62	0.62	0.55	0.60	0.64	0.62	0.66	1.00	0.76
x 17	0.30	0.33	0.31	0.09	0.04	0.09	0.20	0.40	0.42	0.43	0.34	0.39	0.46	0.42	0.50	0.76	1.00

Correlation prediction limit value r ^*get 0.8, known x ₁₆and x ₁₇more weak with the correlativity of other each initializaing variables, but between other initializaing variables, there is stronger correlativity, there is multicollinearity between initializaing variable, need to adopt stepwise regression analysis method to eliminate collinearity.

Forward stepwire regression analytic approach is adopted to screen initializaing variable.Set up the regression model of dependent variable relative to each initializaing variable respectively, calculate the t test statistics of estimates of parameters in each regression model and corresponding probable value p, arrange each initializaing variable by the order of t test statistics (or probable value p from small to large) from big to small.The initializaing variable adding t test statistics large (or probable value p is little) successively carries out successive Regression.Insolation level α gets 5%, when probable value p is less than 5%, retains this initializaing variable, otherwise deletes.Repeat said process, until all initializaing variables are added regression model and completes t inspection.Finally obtain 2 explanatory variables: urban house building completed floor space (is designated as z ₁=x ₁₇) and heat capacity (be designated as z ₂=x ₁).

Step 4, build GM (1,1) based on gray system theory and GM (1,3) model, solving model parameter.

First pre-service is carried out to statistics, calculate the one-accumulate formation sequence of dependent variable and two explanatory variables with

$x_0^{\left(1\right)}=\left(\mathrm{5075,10655,17164,24385,32193,40622,49867,59808}\right)$

$z_1^{\left(1\right)}=\left(\mathrm{1042.2,1778.2,2997,3664.3,4362.5,5297.6,6834.6,8190.7}\right)$

$z_2^{\left(1\right)}=\left(\mathrm{6320,13421,21274,30485.3,40380.6,50909.9,62102.2,73776.5}\right)$

Then to explanatory variable with set up GM (1,1) model respectively:

$\left\{\begin{array}{c}\frac{{\mathrm{dz}}_1^{\left(1\right)}}{\mathrm{dt}}={-a}_{11}{z}_1^{\left(1\right)}+c_1\\ \frac{{\mathrm{dz}}_2^{\left(1\right)}}{\mathrm{dt}}={-a}_{22}{z}_2^{\left(1\right)}+c_2\end{array}\right.$

According to least square method, solve above-mentioned two GM (1,1) model parameter:

${\hat{A}}_1=\left[\begin{array}{c}{\hat{a}}_{11}\\ {\hat{c}}_1\end{array}\right]=\left[\begin{array}{c}0.0139\\ 6.4530\end{array}\right],{\hat{A}}_2=\left[\begin{array}{c}{\hat{a}}_{22}\\ {\hat{c}}_2\end{array}\right]=\left[\begin{array}{c}0.0091\\ 8.7921\end{array}\right]$

Again to dependent variable set up about explanatory variable with gM (1,3) model:

$\frac{{\mathrm{dx}}_0^{\left(1\right)}}{\mathrm{dt}}={-a}_{00}{x}_0^{\left(1\right)}+a_{01}{z}_1^{\left(1\right)}+a_{02}{z}_2^{\left(1\right)}$

According to least square method, solve GM (1,3) model parameter:

${\hat{A}}_0=\left[\begin{array}{c}{\hat{a}}_{00}\\ {\hat{a}}_{01}\\ {\hat{a}}_{02}\end{array}\right]=\left[\begin{array}{c}1.9495\\ 0.2024\\ 1.7536\end{array}\right]$

Step 5, set up gray system state equation, the year heating load of prediction project period.

GM (1, the 1) model of simultaneous two explanatory variables and GM (1, the 3) model of dependent variable, substitute into the model parameter solved in step 4, builds gray system state equation as follows:

$\left[\begin{array}{c}\frac{{\mathrm{dx}}_0^{\left(1\right)}}{\mathrm{d\τ}}\\ \frac{{\mathrm{dz}}_1^{\left(1\right)}}{\mathrm{d\τ}}\\ \frac{{\mathrm{dz}}_2^{\left(1\right)}}{\mathrm{d\τ}}\end{array}\right]=\left[\begin{array}{ccc}-1.9495& 0.2024& 1.7536\\ 0& 0.0139& 0\\ 0& 0& 0.0091\end{array}\right]\left[\begin{array}{c}{x}_0^{\left(1\right)}\\ z_1^{\left(1\right)}\\ z_2^{\left(1\right)}\end{array}\right]+\left[\begin{array}{c}0\\ 6.4530\\ 8.7921\end{array}\right]$

Application Runge-Kutta method solving state equation.If project period is 2011 ~ the year two thousand twenty, by regressive reducing condition non trivial solution, obtain each annual heating load project period:

${\hat{x}}_0^{\left(0\right)}=\left(\begin{array}{c}\mathrm{10359.2,11433.8,12407.8,13591.3,14957.1},\\ \mathrm{16223.9,18005.6,19555.4,21640.3,23735.5}\end{array}\right)$

Harbin City 2003-2012 heating load statistical value and 2011-2020 heating load predicted value are as shown in Figure 2.

Adopt the error ε evaluation and foreca precision of predicted value and statistical value:

$\mathrm{\ϵ}=\frac{|\hat{x}-x|}{x}\×100\%$

Wherein, x is respectively predicted value and statistical value.

The predicted value of Harbin City's heating system year heating load and the error of statistical value of 2011 and 2012 are respectively 1.82%, 8.24%, and be less than acceptable error limit value 10%, precision of prediction can accept.

Claims (8)

1., based on a urban heating system year heating load Forecasting Methodology for gray system theory, it is characterized in that said method comprising the steps of:

Step one, selection initializaing variable, the numerical value over the years of dependent variable: statistics cities and towns year over the years heating load, as dependent variable, adds up town planning index of correlation over the years as initializaing variable;

The degree of association of step 2, the selected initializaing variable of inspection and dependent variable: apply the grey relational grade that grey Relational Analysis Method calculates each initializaing variable and dependent variable, checked the degree of association of selected initializaing variable and dependent variable by the grey relational grade calculated, preliminary screening is carried out to initializaing variable simultaneously;

Step 3, initializaing variable collinearity are checked, and filter out explanatory variable: adopt correlation coefficient matrix inspection initializaing variable collinearity, then by stepwise regression analysis, screening draws explanatory variable;

Step 4, build GM (1,1) based on gray system theory and GM (1, N+1) model, solving model parameter;

Step 5, set up gray system state equation, the year heating load of prediction project period.

2. the urban heating system year heating load Forecasting Methodology based on gray system theory according to claim 1, is characterized in that the computation process of the grey relational grade described in step 2 is as follows:

Step 2 one, employing averaging method carry out nondimensionalization to dependent variable and each initializaing variable:

$X_{i,0}=\frac{x_{i,0}}{{\stackrel{\‾}{x}}_0},X_{i,j}=\frac{x_{i,j}}{{\stackrel{\‾}{x}}_j},(i=\mathrm{1,2},...,m;j=\mathrm{1,2},...,n)---\left(1\right)$

Wherein, X _{i, 0}, x _{i, 0}be respectively dimensionless number and the original value of the i-th annual dependent variable, X _i,j, x _i,jbe respectively dimensionless number and the original value of an i-th annual jth initializaing variable, n is initializaing variable number, and m is the statistics year number of degrees, the average in dependent variable and a jth initializaing variable m year respectively;

Step 2 two, with dimensionless dependent variable for auxiliary sequence X ₀=(X _1,0, X _2,0..., X _{m, 0}) ^t, with dimensionless initializaing variable for subsequence X _j=(X _{1, j}, X _{2, j}..., X _m,j) ^t, calculate each annual auxiliary sequence X ₀with each subsequence X _jabsolute value difference δ _i,j, and get maximal value δ wherein _maxwith minimum value δ _min:

δ _i,j＝|X _i,0-X _i,j| (2)

δ _max＝max{δ _i,j；i＝1,2,…,m；j＝1,2,…,n} (3)

δ _min＝min{δ _i,j；i＝1,2,…,m；j＝1,2,…,n} (4)

Step 2 three, calculate each annual auxiliary sequence X ₀with each subsequence X _jgrey incidence coefficient l _i,j:

$l_{i,j}=\frac{{\mathrm{\δ}}_{\max }+\mathrm{\ζ}\×{\mathrm{\δ}}_{\min }}{{\mathrm{\δ}}_{i,j}+\mathrm{\ζ}\×{\mathrm{\δ}}_{\max }}---\left(5\right)$

Wherein, ζ is resolution ratio;

The grey relational grade of step 2 four, each initializaing variable and dependent variable is:

$L_j=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{l}_{i,j}---\left(6\right).$

3. the urban heating system year heating load Forecasting Methodology based on gray system theory according to claim 2, it is characterized in that described in step 2 to carry out the process of preliminary screening to initializaing variable as follows:

As grey relational grade L _jbe greater than grey relational grade and judge limit value L ^*time, judge that the degree of association of a jth initializaing variable and dependent variable is high, otherwise reject this initializaing variable, complete preliminary screening, wherein L ^*value is 0.6 ~ 0.8.

4. the urban heating system year heating load Forecasting Methodology based on gray system theory according to claim 3, is characterized in that the detailed process of the initializaing variable collinearity inspection described in step 3 is as follows:

Calculate the relative coefficient between initializaing variable:

$r_{\mathrm{uv}}=\frac{\mathrm{Cov}(X_u,X_v)}{\sqrt{\mathrm{Var}\left(X_u\right)}\·\sqrt{\mathrm{Var}\left(X_v\right)}},(u=\mathrm{1,2},...,n;v=\mathrm{1,2},...,n)---\left(7\right)$

Wherein, r _uvbe the relative coefficient between u and v initializaing variable, Var (X _u), Var (X _v) be respectively subsequence X _u, X _vvariance, Cov (X _u, X _v) be subsequence X _u, X _vcovariance;

Relative coefficient between all initializaing variables forms symmetrical correlation coefficient matrix r, and diagonal entry value is 1, more than matrix r diagonal line or the element absolute value of certain row or column lower be designated as | r _pq|, when it is greater than correlation prediction limit value r ^*time, judge u initializaing variable and v initializaing variable correlativity strong, when having strong correlation between at least one initializaing variable and other initializaing variables, judge to there is multicollinearity between initializaing variable, then adopt the collinearity between Stepwise Regression Method elimination initializaing variable, filter out explanatory variable, otherwise all initializaing variables are explanatory variable simultaneously.

5. the urban heating system year heating load Forecasting Methodology based on gray system theory according to claim 4, is characterized in that the detailed process filtering out explanatory variable described in step 3 is as follows:

Stepwise Regression Method is adopted to screen initializaing variable: the horizontal α of specification test, if the probable value p that in regression model, the t test statistics of estimates of parameters is corresponding is less than insolation level α, then retain this initializaing variable, otherwise delete this initializaing variable, the t inspection of all initializaing variables is completed according to forward direction or backward stepwise regression analysis order, the initializaing variable finally retained is the explanatory variable of screening and obtaining, and is designated as z ₁, z ₂..., z _n, wherein α value is 5% ~ 10%.

6. the urban heating system year heating load Forecasting Methodology based on gray system theory according to claim 5, it is characterized in that GM (1, the 1) model of the structure described in step 4 based on gray system theory, the detailed process of solving model parameter is:

The original series of the calendar year statistics value of dependent variable and each explanatory variable composition is designated as with wherein k=1,2 ..., N, carries out the one-accumulate formation sequence obtaining dependent variable and each explanatory variable after adding up, is designated as with

$x_0^{\left(i\right)}=(x_{\mathrm{1,0}},\underset{s=1}{\overset{2}{\mathrm{\Σ}}}{x}_{s,0},...,\underset{s=1}{\overset{i}{\mathrm{\Σ}}}{x}_{s,0},...\underset{s=1}{\overset{m}{\mathrm{\Σ}}}{x}_{s,0})---\left(8\right)$

$z_k^{\left(1\right)}=(z_{1,k},\underset{s=1}{\overset{2}{\mathrm{\Σ}}}{z}_{s,k},...,\underset{s=1}{\overset{i}{\mathrm{\Σ}}}{z}_{s,k},...,\underset{s=1}{\overset{m}{\mathrm{\Σ}}}{z}_{s,k})---\left(9\right)$

Then right set up GM (1,1) model respectively, obtain GM (1, the 1) model of each explanatory variable:

$\frac{{\mathrm{dz}}_k^{\left(1\right)}}{\mathrm{d\τ}}=-a_{\mathrm{kk}}{z}_k^{\left(1\right)}+c_k---\left(10\right)$

Wherein, a _kk, c _kfor the parameter of a kth explanatory variable GM (1,1) model, k=1,2 ..., N, τ are the time;

Discrete differential equation, obtaining matrix equation is:

Y _k＝B _kA _k(11)

Wherein, $Y_k=\left[\begin{array}{c}{z}_{2,k}^{\left(0\right)}\\ z_{3,k}^{\left(0\right)}\\ .\\ .\\ .\\ z_{m,k}^{\left(0\right)}\end{array}\right],$ $B_k=\left[\begin{array}{cc}-\frac{1}{2}[z_{1,k}^{\left(1\right)}+z_{2,k}^{\left(1\right)}]& 1\\ -\frac{1}{2}[z_{2,k}^{\left(1\right)}+z_{3,k}^{\left(1\right)}]& 1\\ .& .\\ .& .\\ .& .\\ -\frac{1}{2}[z_{(m-1),k}^{\left(1\right)}+z_{m,k}^{\left(1\right)}& 1\end{array}\right],$ $A_k=\left[\begin{array}{c}{a}_{\mathrm{kk}}\\ c_k\end{array}\right];$

GM (1, the 1) model parameter of a kth explanatory variable is solved according to least square method:

${\hat{A}}_k=\left[\begin{array}{c}{\hat{a}}_{\mathrm{kk}}\\ {\hat{c}}_k\end{array}\right]={\left(B_k^T{B}_k\right)}^{-1}{B}_k^T{Y}_k---\left(12\right)$

Wherein, for A _k, a _kk, c _kestimated value, k=1,2 ..., N.

7. the urban heating system year heating load Forecasting Methodology based on gray system theory according to claim 6, it is characterized in that the GM (1 of the structure described in step 4 based on gray system theory, N+1) model, the detailed process of solving model parameter is:

To the one-accumulate formation sequence of dependent variable set up GM (1, the N+1) model about N number of explanatory variable one-accumulate formation sequence, obtain GM (1, the N+1) model of dependent variable:

The differential equation form of GM (1, N+1) model is:

$\frac{{\mathrm{dx}}_0^{\left(1\right)}}{\mathrm{d\τ}}={-a}_{00}{x}_0^{\left(1\right)}+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{a}_{0k}{z}_k^{\left(1\right)}---\left(13\right)$

Wherein, a ₀₀, a _0kfor the parameter of GM (1, N+1) model, k=1,2 ..., N;

Discrete differential equation, obtaining matrix equation is:

Y ₀＝B ₀A ₀(14)

Wherein,

$Y_0=\left[\begin{array}{c}{x}_{\mathrm{2,0}}^{\left(0\right)}\\ x_{\mathrm{3,0}}^{\left(0\right)}\\ .\\ .\\ .\\ x_{m,0}^{\left(0\right)}\end{array}\right],$ $B_0=\left[\begin{array}{cccc}-\frac{1}{2}[x_{\mathrm{1,0}}^{\left(1\right)}+x_{\mathrm{2,0}}^{\left(1\right)}]& z_{\mathrm{2,1}}^{\left(1\right)}& ...& z_{2,N}^{\left(1\right)}\\ -\frac{1}{2}[x_{\mathrm{2,0}}^{\left(1\right)}+x_{\mathrm{3,0}}^{\left(1\right)}]& x_{\mathrm{3,1}}^{\left(1\right)}& ...& z_{3,N}^{\left(1\right)}\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ -\frac{1}{2}[x_{(m-1),0}^{\left(1\right)}+x_{m,0}^{\left(1\right)}]& z_{m,1}^{\left(1\right)}& ...& z_{m,N}^{\left(1\right)}\end{array}\right],$ $A_0=\left[\begin{array}{c}{a}_{00}\\ a_{01}\\ .\\ .\\ .\\ a_{0N}\end{array}\right];$

GM (1, N+1) model parameter is solved according to least square method:

${\hat{A}}_0=\left[\begin{array}{c}{\hat{a}}_{00}\\ {\hat{a}}_{01}\\ .\\ .\\ .\\ {\hat{a}}_{0N}\end{array}\right]={\left(B_0^T{B}_0\right)}^{-1}{B}_0^T{Y}_0---\left(15\right)$

Wherein, for A ₀, a ₀₀, a _0kestimated value, k=1,2 ..., N.

8. the urban heating system year heating load Forecasting Methodology based on gray system theory according to claim 7, is characterized in that setting up gray system state equation described in step 5, prediction project period year heating load detailed process be:

The GM (1 of the N number of explanatory variable of simultaneous, 1) model, the GM (1 of dependent variable, N+1) model, formula (12) and formula (15) are substituted into respectively in formula (10) and formula (13), the state equation obtaining gray system is:

The calendar year statistics value of dependent variable and explanatory variable is initial value, the setting year project period number of degrees, and application Runge-Kutta method solving state equation obtains the one-accumulate formation sequence of dependent variable carry out the predicted value that regressive reduction obtains dependent variable again

${\hat{x}}_0=({\hat{x}}_{\mathrm{1,0}}^{\left(1\right)},{\hat{x}}_{\mathrm{2,0}}^{\left(1\right)}-{\hat{x}}_{\mathrm{1,0}}^{\left(1\right)},...,{\hat{x}}_{i,0}^{\left(1\right)}-\underset{s=1}{\overset{i-1}{\mathrm{\Σ}}}{\hat{s}}_{s,0}^{\left(1\right)},...,{\hat{x}}_{M,0}^{\left(1\right)}-\underset{s=1}{\overset{M-1}{\mathrm{\Σ}}}{\hat{x}}_{s,0}^{\left(1\right)}---\left(17\right)$

Wherein, M is the total year number of degrees of statistics and project period.