CN109146196B - 一种居民社区时用水量预测方法 - Google Patents

Abstract

本发明公开了一种居民社区时用水量预测方法，即蜻蜓—高斯过程回归耦合算法，将高斯过程回归与蜻蜓算法进行有机结合，利用擅长全局搜索和对初始值没有依赖的蜻蜓算法进行参数寻优，利用高斯过程回归寻找时用水量训练数据之间的关系，对下一小时的用水量进行预测,给出其分布状况,对居民社区时用水量进行实时预测和区间预测。

Description

一种居民社区时用水量预测方法

技术领域

本发明涉及一种基于改进的蜻蜓—高斯过程回归耦合的居民社区时用水量动态实时区间预测方法，属于水资源管理技术领域。

背景技术

近年来，随着人口的增长、科技的进步、生活水平的提高，水资源诸多问题日渐凸显，水资源日益短缺，对区域水资源进行合理规划、利用和调度势在必行。居民社区用水是水资源消耗的重要组成部分。居民供水系统的优化调度十分重要，不仅要保证用户在不同时间对水量和水压的要求，还要使供水管网压力分布合理，节省大量的泵站电耗，有着十分明显的经济和社会效益。一般来说，供水系统许多运行工作都是以一天(24小时)为周期，其运行方案须至少提前24小时就制定好。而制定该最优运行方案的前提条件就是要准确预测出整个方案期间的时用水量。因此对居民社区时用水量的准确预测至关重要。

针对这一问题，现有的技术方案是：戢钢等人在文章《城市小时级需水量的改进型引力搜索算法-最小二乘支持向量机模型预测》中，用引力搜索算法优化的最小二乘支持向量机模型对城市时用水量进行预测；刘洪波等人在文章《基于人工鱼群神经网络的城市时用水量预测方法》中，利用人工鱼群算法优化的BP神经网络对城市时用水量进行预测；班福忱等人在文章《基于自适应过滤与BP神经网络的城市时用水量组合预测模型》中，通过建立自适应过滤与BP神经网络的组合模型对城市时用水量进行预测；郭泽宇等人在文章《城市用水量组合预测模型及其应用》中，建立了具有季节性时间序列模型和BP神经网络模型优点的组合模型，对城市时用水量进行预测。

经分析，发现现有技术存在如下缺点：1)对于时用水量预测，这些方法都只能得到确定的点预测结果，而居民实际用水量蕴含各种不确定性因素，点预测结果必然存在不同程度的误差，进而使调度决策工作面临一定程度的风险；2)这些方法预测的时用水量是未来某一时刻的值，不能进行实时预测，调度人员无法根据准确的时用水量值进行实时调度。

发明内容

基于此，本发明提出一种新的居民社区时用水量预测方法，即蜻蜓—高斯过程回归耦合算法(DA-GPR，Dragonfly Algorithm-Gaussian Processes Regression)，将高斯过程回归与蜻蜓算法进行有机结合，利用擅长全局搜索和对初始值没有依赖的蜻蜓算法进行参数寻优，利用高斯过程回归(GPR)寻找时用水量训练数据之间的关系，进而对下一小时的用水量进行预测,给出其分布状况,对居民社区时用水量进行实时预测和区间预测。

本发明采用下述的技术方案：

一种基于改进的蜻蜓—高斯过程回归耦合的居民社区时用水量动态实时区间预测方法，包括如下步骤：(1)读取样本数据，产生训练集和测试集；(2)初始参数设置，包括设置蜻蜓种群规模、最大迭代次数和参数l，σ_f，σ_n的取值范围；(3)初始化蜻蜓；(4)将蜻蜓个体位置H的信息依次赋值给参数l，σ_f，σ_n；(5)创建GPR模型，对训练样本进行训练，计算所有蜻蜓个体的目标函数值，并将其作为蜻蜓算法的适应度函数值，计算并保存最优适应度值及其所对应的参数值，种群中的所有蜻蜓个体全部完成分离、对齐、内聚、食物吸引和天敌排斥五种行为操作后，迭代次数加1；(6)判断是否满足终止条件，如果满足条件，输出最优超参数并创建最优的GPR模型；如果不满足条件，则依次更新邻域半径、更新个体位置，然后返回步骤(4)继续迭代；(7)将测试集数据输入创建的最优GPR模型，输出预测值的均值和方差；(8)根据给定置信度，得到时用水量实时区间预测结果。

进一步地，所述初始化蜻蜓包括随机初始化蜻蜓的位置H、步长△H、邻域半径R，惯性权重ω，分离权重s，对齐权重a，内聚权重c，食物吸引权重f，天敌排斥权重e；将蜻蜓个体位置H的信息依次赋值给参数l，σ_f，σ_n是指位置H矩阵的第一行存放参数l的值，第二行存放参数σ_f的值，第三行存放参数σ_n的值，每个蜻蜓个体对应一组参数值。

所述步骤(5)进一步包括：创建GPR模型，对训练样本进行训练，计算所有蜻蜓个体的目标函数值，并将其作为蜻蜓算法的适应度函数值，找出当前的最优个体和最差个体，并将最优个体视为食物，将最差个体视为天敌；蜻蜓每进行一次行为操作，便计算一次当前蜻蜓的最大适应度值：若当前蜻蜓的适应度值大于已保存的最优适应度值，则用当前的蜻蜓的适应度值替代原保存的最优适应度值，作为当前蜻蜓的最优值，并保存当前最优值所对应的蜻蜓的超参数为最优超参数值，否则仍保存原最优适应度值及其所对应的参数值。

优选地，所述终止条件为达到最大迭代次数。

进一步地，目标函数值为

F(x)＝3*MAPE+2*MAE+MSE

其中，MAPE为平均相对误差，MAE为平均绝对误差，MSE为均方误差。

所述超参数为整个高斯过程回归参数集合

其中

是信号方差，M＝diag(l²)，l为方差尺度，σ_n ²为噪声方差。

所述更新邻域半径为邻域半径随迭代次数的增加而成比例地增长；惯性权重ω、分离权重s、对齐权重a、内聚权重c、食物吸引权重f、天敌排斥权重e均为自适应线性递减权重，其最小值为0.4，最大值为0.9；种群规模为10，最大迭代次数均设为20，参数l的取值范围为[-1，1]，参数σ_f的取值范围为[-1，1]，参数σ_n的取值范围为[-1，1]。

区间预测结果表示为：

[m(f_*)-δcov(f_*)，m(f_*)+δcov(f_*)]

其中，f_*为预测值，δ为给定置信水平下预测样本服从后验概率分布的分位数；m(f_*)为均值函数，cov(f_*)为预测值的方差，

其中，K(X,X)＝K_n＝(k_ij)为n×n阶对称正定协方差矩阵；矩阵元素k_ij＝k(x_i,x_j)用来度量x_i和x_j的相关性，x_*为测试点，

为测试点x_*与训练集X之间的协方差阵，k(x_*,x_*)为测试点x_*自身的协方差，I_n是n维单位矩阵。

附图说明

图1是DA-GPR模型区间预测流程图；

图2示出了单日时用水量变化趋势；

图3示出了高斯过程回归(GPR)时用水量预测结果；

图4示出了蜻蜓—高斯过程回归(DA-GPR)时用水量预测结果；

图5示出了改进的蜻蜓—高斯过程回归(DA-GPR)时用水量预测结果；以及

图6示出了不同预测模型预测的相对误差。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

1原理说明：

1.1DA-GPR基本原理

从函数空间角度来看，高斯过程回归的性质由均值函数m(x)和协方差函数k(x,x')确定，定义为：

式中，函数f(x)的高斯过程(Gaussian process,GP)数学表达式为f(x)～GP(m(x),k(x,x')),x,x'∈R^d是任意的随机变量。

对于时用水量回归预测问题，给出如下模型：

y＝f(x)+ε (2)

时用水量实测值y受噪声污染，x是输入向量，

为噪声，则有实测值y和预测值f_*的联合先验分布为：

式中，K(X,X)＝K_n＝(k_ij)为n×n阶对称正定协方差矩阵。矩阵元素k_ij＝k(x_i,x_j)用来度量x_i和x_j的相关性，x_*为测试点，

为测试点x_*与训练集X之间的协方差阵，k(x_*,x_*)为测试点x_*自身的协方差，I_n是n维单位矩阵，σ_n ²为噪声方差。

预测值的后验分布为：

式中，

为预测值的均值，cov(f_*)为预测值的方差，分别为：

本文的协方差函数采用平方指数协方差函数，为：

式中，

是信号方差，M＝diag(l²)，l为方差尺度。整个高斯过程回归参数集合

为超参数。接着用蜻蜓算法对超参数

进行寻优，得到最优超参数。然后联合公式(5)和(6)即可求得测试点x_*所对应的期望和方差。最后，可得到给定置信度下的置信区间作为区间预测结果，可表示为：

[m(f_*)-δcov(f_*)，m(f_*)+δcov(f_*)]。式中，δ为给定置信水平下预测样本服从后验概率分布的分位数。

1.2 DA-GPR目标函数设计

蜻蜓—高斯过程回归耦合时用水量预测模型的最终目标是对未来的时用水量进行较准确的动态实时区间预测，即具有较好的泛化能力。因此在使用蜻蜓算法优化高斯过程回归的参数时，选取能代表泛化能力的目标函数是非常重要的。常规的蜻蜓优化算法将训练集上获得最小平均相对误差(MAPE)的参数作为最优解。即，将训练集上获得最小平均相对误差(MAPE)作为目标函数，对蜻蜓个体进行选择，目标函数值越小，蜻蜓个体越优，对应的参数越优。由此可见，目标函数的选择至关重要。

居民社区每小时的用水量比较少，基本上在几立方米到十几立方米，仅仅采用平均相对误差(MAPE)作为目标函数值进行参数寻优的效果不理想。在此基础上，本文对DA-GPR进行改进，将体现绝对误差大小的平均绝对误差(MAE)和均方误差(MSE)引入到目标函数中，将平均相对误差(MAPE)、平均绝对误差(MAE)和均方误差(MSE)进行线性组合，构造出新的目标函数，如式(8)。

F(x)＝3*MAPE+2*MAE+MSE (8)

1.3 DA-GPR局部搜索公式设计

蜻蜓群体在一个小的区域内来回飞翔寻找猎物以及躲避外敌。通过模拟蜻蜓群体之间寻找猎物以及躲避外敌的社会互动行为来实现优化的过程，即对算法的局部搜索过程。

蜻蜓个体在群体运动中可以分为分离、对齐、内聚、食物吸引和天敌排斥等五种行为模式。该行为的具体意义和数学表达方法如下：

(1)分离：指蜻蜓与相邻个体之间避免碰撞。该行为的计算公式如下:

(2)对齐：指相邻个体之间倾向于保持相同的速度。该行为的计算公式如下:

(3)内聚：指蜻蜓个体试图向自己认为所属的群体靠近。该行为的计算公式如下:

(4)食物吸引力：指食物对蜻蜓的吸引力。该行为的计算公式如下:

F_i＝H⁺-H (12)

(5)天敌排斥力：指蜻蜓对天敌的排斥力。该行为的计算公式如下:

E_i＝H^-+H (13)

式中，S_i表示第i个个体的分离量；H是当前个体的位置；H_j是相邻个体j的位置；K为相邻个体的数量；A_i表示第i个个体的对齐量；V_j为第j个相邻个体的速度；C_i表示第i个个体的内聚量；F_i表示第i个个体对食物的吸引力；H⁺表示食物所在的位置；E_i表示第i个个体需逃离天敌的距离；H^-表示天敌所在的位置。

1.4DA-GPR全局搜索公式设计

在蜻蜓寻优的过程中，大量的蜻蜓群体朝着同一个方向进行长距离迁移，这个群体活动就是对算法的全局搜索过程。

根据5种蜻蜓行为，下一代蜻蜓的步长计算如下：

ΔH_t+1＝(sS_i+aA_i+cC_i+fF_i+eE_i)+ωΔH_t (14)

有邻近蜻蜓时，下一代蜻蜓的位置为：

H_t+1＝H_t+ΔH_t+1 (15)

无邻近蜻蜓时，设置为随机游走(Le'vy飞行)行为，下一代蜻蜓的个体位置为:

H_t+1＝H_t+Le'vy(d)×H_t (16)

式(14)(15)和(16)中，t表示当前迭代次数；i表示第i个蜻蜓个体；H_t表示当前第t代种群个体位置；ΔH_t+1表示第t+1代种群位置更新步长；H_t+1表示第t+1代种群个体位置；s，a，c，f，e分别对应于5种行为的权重；ω是惯性权重；d是个体位置向量的维数。

Le'vy函数计算公式如下：

式(17)和(18)中，r₁，r₂是[0,1]随机数，Γ(w)＝(w-1)！，β是一个常数(这里取为0.5)。蜻蜓算法在寻优的过程中，对每个个体的邻近个体数量的计算是非常重要的，因此这里假定一个邻域半径，该半径随迭代次数的增加而成比例地增长。同时，为了达到算法的局部搜索和全局搜索的平衡，权重(s，a，c，f，e和_ω)在优化过程中自适应地调整。

1.5DA-GPR算法步骤

流程图如图1所示。具体实现步骤如下：

1.读取样本数据，产生训练集和测试集。

2.初始参数设置：设置蜻蜓种群规模、最大迭代次数和参数l，σ_f，σ_n的取值范围。

3.初始化蜻蜓：随机初始化蜻蜓的位置H、步长△H、邻域半径R，惯性权重ω，分离权重s，对齐权重a，内聚权重c，食物吸引权重f，天敌排斥权重e。

4.将蜻蜓个体位置H的信息依次赋值给l，σ_f，σ_n，其中位置H矩阵的第一行存放参数l的值，第二行存放参数σ_f的值，第三行存放参数σ_n的值，每个蜻蜓个体对应一组参数值。

5.创建GPR模型，对训练样本进行训练，计算所有蜻蜓个体的目标函数值，并将其作为蜻蜓算法的适应度函数值，找出当前的最优个体和最差个体，并将最优个体视为食物，将最差个体视为天敌。蜻蜓每进行一次行为操作，便计算一次当前蜻蜓的最大适应度值：若当前蜻蜓的适应度值大于已保存的最优适应度值，则用当前的蜻蜓的适应度值替代原保存的最优适应度值，作为当前蜻蜓的最优值，并保存当前最优值所对应的蜻蜓的超参数为最优超参数值，否则仍保存原最优适应度值及其所对应的参数值。种群中的所有蜻蜓个体全部完成分离、对齐、内聚、食物吸引和天敌排斥等五种行为操作后，迭代次数加1。

6.判断是否满足终止条件，如果满足条件，输出最优超参数并创建最优的GPR模型；如果不满足条件，则依次更新邻域半径、更新个体位置，然后返回步骤4继续迭代。

7.将测试集数据输入创建的最优GPR模型，输出预测值的均值和方差。

8.根据给定置信度，得到时用水量实时区间预测结果。

步骤6中的终止条件优选地可以是达到最大迭代次数。

2仿真实验

2.1仿真实验平台

本次实验的运行环境为：操作系统：Windows 7；处理器：Intel(R)Core(TM)i7-2640M；主频：2.80GHz；内存：4.00GB；编程工具：MATLAB R2015(b)。

2.2样本数据

本文所用的数据来自河北工程大学在线水量检测平台，采用居民家属院2016年6月17日到7月6日，共21天每天24小时的实际社区时用水量资料，如图2所示。将样本集21天每小时用水量中的前20天480组时用水量资料作为训练集，最后一天24小时用水量资料作为测试集。

2.3模型输入和输出

从图2中可以明显看出居民社区时用水量具有明显的24小时周期性，为预测下一时刻用水量，模型的输入取该时刻前24小时的时用水量。也就是说，模型输出是居民社区下一时刻用水量，模型输入是该时刻前24小时的时用水量。

2.4参数设置

惯性权重ω、分离权重s、对齐权重a、内聚权重c、食物吸引权重f、天敌排斥权重e均为自适应线性递减权重，其最小值为0.4，最大值为0.9；种群规模为10，最大迭代次数均设为20，参数l的取值范围为[-1，1]，参数σ_f的取值范围为[-1，1]，参数σ_n的取值范围为[-1，1]。

2.5对比模型及评价标准

为了使DA-GPR模型的预测结果具有可比性，选择目标函数改进前的DA-GPR模型和共轭梯度法选择参数的GPR模型作为对比模型，模型性能采用平均相对误差(MAPE)、平均绝对误差(MAE)和均方误差(MSE)进行衡量，它们的定义为：

式中：N为测试集样本的大小，y_i为实测值，y_i为预测值。

为了验证本文创建的改进的蜻蜓—高斯过程回归耦合模型的预测可靠性，本文将其与高斯过程回归和蜻蜓—高斯过程回归两种预测模型进行试验对比。在95％置信度下，不同模型的预测结果如图3、图4和图5。从图3可知，高斯过程回归预测模型的预测误差较大，极个别点落在了预测区间外。从图4可知，蜻蜓—高斯过程回归预测模型预测的结果走势基本上与实测值走势一致但是部分点的预测误差较大。从图5可知，改进的蜻蜓—高斯过程回归预测模型预测的结果走势与实测值一致，预测值与实测值基本上重合，误差较小。

不同预测模型的预测相对误差如图6。从图6可以看出，改进的蜻蜓—高斯过程回归预测模型相对误差最大为0.019，蜻蜓—高斯过程回归预测模型相对误差最大为0.145，,高斯过程回归预测模型相对误差最大为1.653。显然，改进的蜻蜓—高斯过程回归预测模型具有更高的预测精度。

采用平均相对误差(MAPE)、平均绝对误差(MAE)和均方误差(MSE)对不同模型的预测结果进行分析，结果如表1所示。通过表1可以看出，改进的蜻蜓—高斯过程回归模型的居民社区时用水量预测精度优于蜻蜓—高斯过程回归模型和高斯过程回归模型，取得了较好的预测效果。

表1不同预测模型误差分析

本申请将高斯过程回归与蜻蜓算法进行有机结合，提出一种新的居民社区时用水量预测方法，即蜻蜓—高斯过程回归耦合算法(DA-GPR)，对居民社区时用水量进行预测；本申请所述的蜻蜓—高斯过程回归耦合预测方法与常规方法相比，不仅能够预测未知量的期望值，还能给出其分布状况即预测波动区间,同时也能进行实时动态预测。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于改进的蜻蜓—高斯过程回归耦合(Dragonfly Algorithm-GaussianProcesses Regression,DA-GPR)的居民社区时用水量动态实时区间预测方法，包括如下步骤：

(1)读取样本数据，产生训练集和测试集；

(2)初始参数设置，设置蜻蜓种群规模、最大迭代次数和参数l，σ_f，σ_n的取值范围；

(3)初始化蜻蜓；

(4)将蜻蜓个体位置H的信息依次赋值给参数l，σ_f，σ_n，其中l为方差尺度，

是信号方差，σ_n ²为噪声方差；

(5)创建高斯过程回归(Gaussian Processes Regression,GPR)模型，对训练样本进行训练，计算所有蜻蜓个体的目标函数值，并将其作为蜻蜓算法的适应度函数值，计算并保存最优适应度值及其所对应的参数值，种群中的所有蜻蜓个体全部完成分离、对齐、内聚、食物吸引和天敌排斥五种行为操作后，迭代次数加1；

(6)判断是否满足终止条件，如果满足条件，输出最优超参数并创建最优的GPR模型；如果不满足条件，则依次更新邻域半径、更新个体位置，然后返回步骤(4)继续迭代；

(7)将测试集数据输入创建的最优GPR模型，输出预测值的均值和方差；

(8)根据给定置信度，得到时用水量实时区间预测结果。

2.根据权利要求1所述的方法，所述初始化蜻蜓包括随机初始化蜻蜓的位置H、步长△H、邻域半径R，惯性权重ω，分离权重s，对齐权重a，内聚权重c，食物吸引权重f，天敌排斥权重e。

3.根据权利要求2所述的方法，所述将蜻蜓个体位置H的信息依次赋值给参数l，σ_f，σ_n是指位置H矩阵的第一行存放参数l的值，第二行存放参数σ_f的值，第三行存放参数σ_n的值，每个蜻蜓个体对应一组参数值。

4.根据权利要求3所述的方法，所述步骤(5)进一步包括：创建GPR模型，对训练样本进行训练，计算所有蜻蜓个体的目标函数值，并将其作为蜻蜓算法的适应度函数值，找出当前的最优个体和最差个体，并将最优个体视为食物，将最差个体视为天敌；蜻蜓每进行一次行为操作，便计算一次当前蜻蜓的最大适应度值：若当前蜻蜓的适应度值大于已保存的最优适应度值，则用当前的蜻蜓的适应度值替代原保存的最优适应度值，作为当前蜻蜓的最优值，并保存当前最优值所对应的蜻蜓的超参数为最优超参数值，否则仍保存原最优适应度值及其所对应的参数值。

5.根据权利要求4所述的方法，所述终止条件为达到最大迭代次数。

6.根据权利要求1或4所述的方法，所述目标函数值为

F(x)＝3*MAPE+2*MAE+MSE

其中，MAPE为平均相对误差，MAE为平均绝对误差，MSE为均方误差。

7.根据权利要求6所述的方法，所述超参数为整个高斯过程回归参数集合

其中M＝diag(l²)。

8.根据权利要求7所述的方法，所述更新邻域半径为邻域半径随迭代次数的增加而成比例地增长；惯性权重ω、分离权重s、对齐权重a、内聚权重c、食物吸引权重f、天敌排斥权重e均为自适应线性递减权重，其最小值为0.4，最大值为0.9；种群规模为10，最大迭代次数均设为20，参数l的取值范围为[-1，1]，参数σ_f的取值范围为[-1，1]，参数σ_n的取值范围为[-1，1]。

9.根据权利要求8所述的方法，所述区间预测结果表示为：[m(f_*)-δcov(f_*)，m(f_*)+δcov(f_*)]

其中，f_*为预测值，δ为给定置信水平下预测样本服从后验概率分布的分位数；m(f_*)为均值函数，cov(f_*)为预测值的方差，

其中，K(X,X)＝K_n＝(k_ij)为n×n阶对称正定协方差矩阵；矩阵元素k_ij＝k(x_i,x_j)用来度量x_i和x_j的相关性，x_*为测试点，

为测试点x_*与训练集X之间的协方差阵，k(x_*,x_*)为测试点x_*自身的协方差，I_n是n维单位矩阵。

Images