基于灰色系统理论的城镇供热系统年供热量预测方法

技术领域

本发明涉及一种城镇供热系统年供热量预测方法，尤其涉及一种基于灰色系统理论的 城镇供热系统年供热量预测方法，属于供热系统负荷预测技术领域。

背景技术

随着我国城镇化进程和人民生活水平的提高，对建筑室内环境的舒适度和可控制性提 出了更高的要求，我国城镇供热面积和供热系统供热量不断增加，1990年至2012年，我 国集中供热面积从2.13亿平方米增长至51.84亿平方米，年供热总量从2.87亿吉焦增加 到29.5亿吉焦，供热系统能耗已成为社会总能耗的重要组成部分。但在我国紧张的能源 供应现状和大力发展低碳经济的背景下，需合理规划用能结构，提高能源效率。

城市能源系统规划建立在对能源需求合理预测的基础上，供热规划主要分为两类：一 类是针对单一供热系统，基于近期和远期扩建规模，预测中、长期供热负荷，规划系统扩 建容量；另一类是针对城镇，在原有分散或集中供热基础上，考虑城镇的近期和远期规划 发展，预测供热负荷或年供热量，规划城镇中全部供热系统的能源产出与供应。目前国内 外针对单一供热系统的负荷预测往往基于对象系统的热力和水力工况的运行参数、室外气 象参数等历史数据，无法用于城镇供热系统负荷预测，而在城镇供热系统规划时，主要采 用指标概算法或度日法预测热负荷或年供热量，估算依据来源单一，不能反映城镇的经济、 人口、产业和建设等发展水平。

因此，亟需一种可应用于城镇能源系统或供热系统规划，并与城镇规划依据一致，不 仅能够体现供热工程发展水平，亦能够反映城镇经济、人口、产业和建设等发展水平的供 热系统年供热量预测方法。

发明内容

本发明的目的是提出一种基于灰色系统理论的城镇供热系统年供热量预测方法，以解 决针对单一供热系统的预测方法不能用于城镇供热系统规划的问题。

本发明为解决上述技术问题所采用的技术方案是：

本发明所述的基于灰色系统理论的城镇供热系统年供热量预测方法，是按照以下步骤 实现的：

步骤一、选择初始变量，因变量的历年数值：统计城镇历年年供热量作为因变量，统 计城镇规划历年的相关指标作为初始变量；

步骤二、检验所选初始变量与因变量的关联度：应用灰色关联度分析法计算各初始变 量与因变量的灰色关联度，通过计算出的灰色关联度检验所选初始变量与因变量的关联 度，同时对初始变量进行初步筛选；

步骤三、初始变量共线性检验，并筛选出解释变量：采用相关性系数矩阵检验初始变 量共线性，再通过逐步回归分析，筛选得出解释变量；

步骤四、构建基于灰色系统理论的GM(1，1)和GM(1，N+1)模型，求解模型参数；

步骤五、建立灰色系统状态方程，预测规划期的年供热量。

本发明的有益效果是：

1.本发明选取的解释变量与供热系统运行参数无关，不受城镇中供热系统运行参数 的限制，而与城镇发展状况密切相关，可反映城镇供热工程、经济、人口、产业和建设等 发展水平，与城镇规划依据一致。

2.本发明以城镇为供热范围，不受城镇中供热系统的数量限制，预测的年供热量反 映整个城镇供热负荷，同时根据气象参数，可将年供热量折算为供暖室外计算温度下的热 负荷，作为城镇热源规划的依据。

3.本发明通过筛选消除初始变量间的相关性，精简为若干解释变量，可简化预测模 型，同时保证对因变量的解释性。

4.本发明基于灰色系统理论建立的状态方程适于样本数据较少的城镇年供热量的预 测，同时本发明具备学习性，可根据逐年的数据更新，重新筛选解释变量并计算模型参数， 建立新的状态方程，保证预测的实时性与可靠性。

附图说明

图1为本发明方法的流程图；

图2为本发明实施例中哈尔滨市2003-2012年供热量统计值与2011-2020年供热量预 测值示意图。

具体实施方式

结合附图进一步详细说明本发明的具体实施方式。

具体实施方式一：下面结合图1说明本实施方式，本实施方式所述的一种基于灰色系 统理论的城镇供热系统年供热量预测方法，其特征在于所述方法包括以下步骤：

步骤一、选择初始变量，因变量的历年数值：以城镇供热系统年供热量为预测对象， 统计城镇历年年供热量作为因变量，由于供热系统作为北方城镇基础设施建设的一部分， 除供热工程自身的发展水平外，还受到城镇经济、人口、产业和建设等发展水平的制约， 因此统计城镇规划历年的相关指标作为初始变量；

所述相关指标为：供热能力、供热面积、供热管道长度、居住用地、公共设施用地、 工业用地、城市年底总人口、城市总户数、地区生产总值、城市人均地区生产总值、城市 第一产业生产总值、城市第二产业生产总值、城市第三产业生产总值、居民家庭人均可支 配收入、城镇居民人均消费支出，城市房屋建筑施工面积和城市房屋建筑竣工面积，共 17个初始变量的统计值。

步骤二、检验所选初始变量与因变量的关联度：应用灰色关联度分析法计算各初始变 量与因变量的灰色关联度，通过计算出的灰色关联度检验所选初始变量与因变量的关联 度，同时对初始变量进行初步筛选，剔除灰色关联度较低的初始变量；

步骤三、初始变量共线性检验，并筛选出解释变量：由于初始变量之间可能存在的共 线性会影响预测模型的准确性。因此采用相关性系数矩阵检验初始变量共线性，再通过逐 步回归分析，在兼顾对因变量解释性和模型精度的前提下，筛选得出解释变量；

步骤四、构建基于灰色系统理论的GM(1，1)和GM(1，N+1)模型，求解模型参数；

步骤五、建立灰色系统状态方程，预测规划期的年供热量。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤二所述的灰色关联度 的计算过程如下：

步骤二一、采用均值法对因变量和各初始变量进行无量纲化：

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

其中，X_i，0、x_i，0分别为第i年度因变量的无量纲值和原始值，X_i，j、x_i，j分别为第i年度 第j个初始变量的无量纲值和原始值，n为初始变量个数，m为统计年度数，分别 是因变量和第j个初始变量m年度的均值；

步骤二二、以无量纲因变量为母序列X₀＝(X_1，0，X_2，0，…，X_m，0)^T，以无量纲初始变量为子 序列X_j＝(X_1，j，X_2，j，…，X_m，j)^T，计算各年度母序列X₀与各子序列X_j的绝对值差δ_i，j，并取其中 的最大值δ_max和最小值δ_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)

步骤二三、计算各年度母序列X₀与各子序列X_j的灰色关联系数l_i，j：

$l_{i,j}=\frac{{\delta }_{\max }+{\zeta \times \delta }_{\min }}{{\delta }_{i,j}+{\zeta \times \delta }_{\max }}---\left(5\right)$

其中，ζ为分辨系数，其值为0.5；

步骤二四、各初始变量与因变量的灰色关联度为：

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

具体实施方式三：本实施方式与具体实施方式一或二不同的是：

步骤二所述的对初始变量进行初步筛选的过程如下：

当灰色关联度L_j大于灰色关联度判定限值L^*时，判定第j个初始变量与因变量的关 联度高，否则剔除该初始变量，完成初步筛选，其中L^*取值为0.6～0.8。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤三所述的初 始变量共线性检验的具体过程如下：

计算初始变量间的相关性系数：

$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)$

其中，r_uv为第u个和第v个初始变量间的相关性系数，Var(X_u)、Var(X_v)分别为子序 列X_u、X_v的方差，Cov(X_u，X_v)为子序列X_u、X_v的协方差；

所有初始变量间的相关性系数构成对称的相关性系数矩阵r，且对角线元素值为1， 矩阵r对角线以上或下某行或列的元素绝对值记为|r_pq|，当其大于相关性判定限值r^*时， 判定第u个初始变量与第v个初始变量相关性强，当至少一个初始变量与其他初始变量之 间有强相关性时，判定初始变量之间存在多重共线性，则采用逐步回归分析法消除初始变 量间的共线性，同时筛选出解释变量，否则所有初始变量即为解释变量，其中r^*取值0.8。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤三所述的筛 选出解释变量的具体过程如下：

采用逐步回归分析法对初始变量进行筛选：设定检验水平α，若回归模型中参数估计 值的t检验统计量对应的概率值p小于检验水平α，则保留该初始变量，否则删除该初始 变量，按照前向或后向逐步回归分析顺序完成所有初始变量的t检验，最后保留的初始变 量即为筛选得到的解释变量，记为z₁，z₂，…，z_N，其中α取值为5％～10％。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤四所述的构 建基于灰色系统理论的GM(1，1)模型，求解模型参数的具体过程为：

首先对统计数据进行预处理，将因变量和各解释变量的历年统计值组成的原始序列记 为和其中k＝1，2，…，N，进行累加后得到因变量和 各解释变量的一次累加生成序列，记为和

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

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

然后对分别建立GM(1，1)模型，得到各解释变量的GM(1，1)模型：

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

其中，a_kk、c_k为第k个解释变量GM(1，1)模型的参数，k＝1，2，…，N，τ为时间；

离散微分方程，得到矩阵方程为：

Y_k＝B_kA_k               (11)

其中，$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);$

根据最小二乘法求解第k个解释变量的GM(1，1)模型参数：

${\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)$

其中，为A_k、a_kk、c_k的估计值，k＝1，2，…，N。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：步骤四所述的构建基 于灰色系统理论的GM(1，N+1)模型，求解模型参数的具体过程为：

对因变量的一次累加生成序列建立关于N个解释变量一次累加生成序列的 GM(1，N+1)模型，得到因变量的GM(1，N+1)模型：

GM(1，N+1)模型反映N个解释变量和因变量本身对于因变量一阶导数的影响。其微分 方程形式为：

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

其中，a₀₀、a_0k为GM(1，N+1)模型的参数，k＝1，2，…，N；

离散微分方程，得到矩阵方程为：

Y₀＝B₀A₀              (14)

其中，

$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)模型参数：

${\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)$

其中，为A₀、a₀₀、a_0k的估计值，k＝1，2，…，N。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：步骤五所述的建立灰 色系统状态方程，预测规划期的年供热量的具体过程为：

联立N个解释变量的GM(1，1)模型(公式(10))，因变量的GM(1，N+1)模型(公式(13))， 将公式(12)和公式(15)分别代入公式(10)和公式(13)中，得到灰色系统的状态方程为：

因变量和解释变量的历年统计值为初值，设定规划期年度数，应用Runge-Kutta法求 解状态方程得到因变量的一次累加生成序列再进行累减还原得到因变量的预测值

${\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}{\Sigma }}{\hat{x}}_{s,0}^{\left(1\right)},...,{\hat{x}}_{M,0}^{\left(1\right)}-\underset{s=1}{\overset{M-1}{\Sigma }}{\hat{x}}_{s,0}^{\left(1\right)})---\left(17\right)$

其中，M为统计与规划期的总年度数。

本发明实施例如下：

以下结合哈尔滨市年供热量预测进一步阐明。

为验证本发明所述预测方法的精度，取2003～2010年统计数据预测2011～2020年的年 供热量，其中2011和2012年供热量预测值与统计值进行对比，计算预测误差。

步骤一、选择初始变量，因变量的历年数值。

根据《中国城乡建设数据库》摘录哈尔滨市2003-2012年供热量(x₀)作为因变量统 计值，见表1。

表1哈尔滨市2003-2012年供热量统计值

根据《中国城市建设统计年鉴》、《中国区域经济统计年鉴》和《中国城乡建设数据库》 摘录哈尔滨市2003-2010供热能力(x₁)、供热面积(x₂)、供热管道长度(x₃)、居住用地 (x₄)、公共设施用地(x₅)、工业用地(x₆)、城市年底总人口(x₇)、城市总户数(x₈)、 地区生产总值(x₉)、城市人均地区生产总值(x₁₀)、城市第一产业生产总值(x₁₁)、城市 第二产业生产总值(x₁₂)、城市第三产业生产总值(x₁₃)、居民家庭人均可支配收入(x₁₄)、 城镇居民人均消费支出(x₁₅)、城市房屋建筑施工面积(x₁₆)和城市房屋建筑竣工面积(x₁₇) 等17个初始变量的统计值，见表2。

表2哈尔滨市2003-2010初始变量统计值

步骤二、检验所选初始变量与因变量的关联度。

采用均值法对因变量x₀和初始变量x₁～x₁₇进行无量纲化，构成母序列X₀与子序列X_j， 计算2003～2010年度各子序列与母序列的绝对值差，得到其最大值δ_max＝0.40862，最小值 δ_min＝0.00146。

分辨系数ζ取0.5，计算2003～2010年度母序列X₀与各子序列X_j的灰色关联系数l_i，j， 得到各初始变量与因变量的灰色关联度L_j，见表3。

表3初始变量x₁～x₁₇与因变量x₀的灰色关联度L₁～L₁₇

灰色关联度判定限值L^*取0.6，各初始变量与因变量的关联度L_j均高于0.6，说明所 选17个初始变量与因变量关联度高。

步骤三、初始变量共线性检验，并筛选出解释变量。

计算初始变量间相关性系数，构成相关性系数矩阵r，见表4。

表4初始变量的相关性系数矩阵r

相关性判定限值r^*取0.8，可知x₁₆和x₁₇与其他各初始变量的相关性较弱，但是其他 初始变量间具有较强的相关性，初始变量之间存在多重共线性，需要采用逐步回归分析方 法消除共线性。

采用前向逐步回归分析法对初始变量进行筛选。分别建立因变量相对于每个初始变量 的一元回归模型，计算各回归模型中参数估计值的t检验统计量和相应的概率值p，按t 检验统计量从大到小(或概率值p从小到大)的顺序排列各初始变量。依次加入t检验统 计量大(或概率值p小)的初始变量进行逐步回归。检验水平α取5％，当概率值p小于 5％时，保留该初始变量，否则删除。重复上述过程，直到把所有初始变量加入回归模型 并完成t检验。最后得到2个解释变量：城市房屋建筑竣工面积(记为z₁＝x₁₇)和供热能 力(记为z₂＝x₁)。

步骤四、构建基于灰色系统理论的GM(1，1)和GM(1，3)模型，求解模型参数。

首先对统计数据进行预处理，计算因变量和两个解释变量的一次累加生成序列和

$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)$

然后对解释变量和分别建立GM(1，1)模型：

$\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)$

根据最小二乘法，求解上述两个GM(1，1)模型参数：

${\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)$

再对因变量建立关于解释变量和的GM(1，3)模型：

$\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)}$

根据最小二乘法，求解GM(1，3)模型参数：

${\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)$

步骤五、建立灰色系统状态方程，预测规划期的年供热量。

联立两个解释变量的GM(1，1)模型和因变量的GM(1，3)模型，代入步骤四中求解的模 型参数，构建灰色系统状态方程如下：

$\left(\begin{array}{c}\frac{{\mathrm{dx}}_0^{\left(1\right)}}{\mathrm{d\tau }}\\ \frac{{\mathrm{dz}}_1^{\left(1\right)}}{\mathrm{d\tau }}\\ \frac{{\mathrm{dz}}_2^{\left(1\right)}}{\mathrm{d\tau }}\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)$

应用Runge-Kutta法求解状态方程。设规划期为2011～2020年，通过累减还原状态方 程的解，得到规划期各年度供热量：

${\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)$

哈尔滨市2003-2012年供热量统计值与2011-2020年供热量预测值如图2所示。

采用预测值与统计值的误差ε评价预测精度：

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

其中，x分别为预测值和统计值。

2011和2012年的哈尔滨市供热系统年供热量的预测值与统计值的误差分别为 1.82％、8.24％，小于可接受误差限值10％，预测精度可接受。