CN111144644B

CN111144644B - 基于变分异方差高斯过程回归的短期风速预测方法

Info

Publication number: CN111144644B
Application number: CN201911348832.1A
Authority: CN
Inventors: 张楚; 彭甜; 夏鑫; 王业琴; 赵环宇; 孙娜; 张涛
Original assignee: Huaiyin Institute of Technology
Current assignee: Huaiyin Institute of Technology
Priority date: 2019-12-24
Filing date: 2019-12-24
Publication date: 2022-01-28
Anticipated expiration: 2039-12-24
Also published as: CN111144644A

Abstract

本发明涉及短期风速预测技术领域，公开了一种基于变分异方差高斯过程回归的短期风速预测方法，首先采用具适应性噪声的完全集合经验模态分解将原始风速时间序列分解为规律性更强的子序列集合，计算每个子序列的偏自相关函数值，并选择95%置信水平下显著的滞后时间序列作为输入变量；然后，采用变分异方差高斯过程回归对每个子序列进行训练与预测，最后，将所有子序列的预测结果进行合并，得到风速时间序列的最终预测结果。与现有技术相比，本发明采用变分异方差高斯过程回归模型对风速时间序列进行预测，其预测能力更强，且变分异方差高斯过程回归模型性能优于标准高斯过程回归模型，能够获得更高的预测精度。

Description

基于变分异方差高斯过程回归的短期风速预测方法

技术领域

本发明涉及短期风速预测领域，特别涉及一种基于变分异方差高斯过程回归的短期风速预测方法。

背景技术

随着世界范围内化石燃料的减少和环境污染水平的恶化，可再生能源越来越受到人们的重视。作为一种清洁、蕴量巨大、可循环利用的可再生能源，风能已成为电力系统最具发展前景的能源之一。但由于风速的间歇性、波动性和随机性，使得风能并网的难度增大。准确可靠的短期风速预测方法不仅能为经济可靠的发电计划制定提供保障，而且能够减轻风电并网难度，提高风力发电能源开发利用效率。

根据不同风速预测模型的机理，风速预测方法主要可分为物理模型、时间序列模型、智能预测模型三类。物理模型主要基于数值天气预报NWP模型。NWP模型通常采用风速、温度、湿度、地形粗糙度和障碍物等来建立一系列的力学和热力学方程，适用于未来三天的中长期风速预报。然而，NWP预测模型建模过程复杂，计算成本高，不适合于短期风速预测。与物理模型不同，时间序列和智能预测模型的建立主要基于历史观测风速数据。时间序列模型通常假设时间序列是平稳的并且服从正态分布。然而，实际风速时间序列一般具有较强的随机性和波动性。因此，在处理强非线性风速时间序列时，时间序列模型很难获得良好的结果。

智能预测模型主要包括神经网络、支持向量机、高斯过程回归等。其中，高斯过程回归是一种新型的基于统计和贝叶斯理论的智能回归模型，它同时具备贝叶斯网络的数学推理能力和支持向量机的自适应能力，尤其适合处理小样本、高维、强非线性的复杂线性回归问题。然而，标准高斯过程回归模型假设噪声符合高斯分布，且整个数据集的噪声方差为常数，使得实际风速的预测误差较大，且具有很强的随机性和波动性。异方差高斯过程回归是标准高斯过程回归的一个扩展，它采用两个高斯过程分别对无噪声输出和噪声进行建模。异方差高斯过程回归能够较好的捕捉风速的随机性和波动性，然而，异方差高斯过程回归模型的最大后验概率估计不能将所有潜在的自变量结合起来，容易遇到过拟合问题。

发明内容

发明目的：针对风电并网时风速的强波动性对电能质量和电力系统安全稳定运行造成的不利影响，本发明提出了基于变分异方差高斯过程回归的短期风速预测方法，得到风速时间序列的最终预测结果，以提高风速时间序列的预测精度。

技术方案：本发明提供了一种基于变分异方差高斯过程回归的短期风速预测方法，包括如下步骤：

步骤一：收集风电场历史实测短期风速数据，根据历史实测风速数据建立风速时间序列，通过完全集合经验模态分解将其分解为K个本证模态函数分量IMK(k＝1,2,3,...,K)和一个残差变量R_k；

步骤二：对于步骤一中每个本征模态函数分量IMF_k(k＝1,2,3,...,K)和残差R_k，计算其偏自相关函数值，选择95％置信水平下显著的滞后时间序列作为输入变量，目标向量为输出变量，并将全部样本数据的前70％作为训练样本，剩余30％作为检验样本；

步骤三：对于步骤一中每个本征模态函数分量IMF_k(k＝1,2,3,...,K)和残差R_k，建立变分异方差高斯过程回归模型，输出最优变分参数Λ和模最优型超参数

步骤四：将步骤三所得的最优变分参数Λ和最优模型超参数

以及步骤二中的检验样本代入变分异方差高斯过程回归预测模型，得到检验阶段的预测值，并对所有子模态的预测结果进行求和重构，得到原始风速时间序列的预测值；

步骤五：计算预测值和实测值的均方根误差RMSE、平均绝对误差MAE、相关系数R和平均绝对百分比误差MAPE，评价所提出的风速预测模型的性能。

进一步地，所述变分异方差高斯过程回归模型是引入边缘变分方法估计异方差高斯过程回归的预测分布，所述异方差高斯过程记为G3，其由两个方差相同的高斯过程组成：第一个高斯过程模型称为G1，对无噪声目标输出进行建模；第二个高斯过程模型称为G2，对数噪声进行建模，即，z(x_i)＝lg(r(x_i))，该z过程由协方差函数k_z控制，其超参数为θ_z。

进一步地，所述步骤三的变分异方差高斯过程回归的步骤如下：

3.1)对于给定的训练集D＝{(x_i,y_i)(i＝1,2,...,n)}，应用第一个高斯过程G1预测输出向量y，构建GPR模型，其中x_i表示d维输入向量，y_i表示对应的输出向量，n表示训练样本的数量；

3.2)根据G1估计训练数据的噪声水平，即z_i＝lg{var[y_i,G(x_i,X,y)]}，形成新的数据集D'＝{(x_i,z_i)(i＝1,2,...,n)}；

3.3)运用第二个高斯过程G2对数据集D′建模；

3.4)运用异方差高斯过程G3对数据集D建模，其中，用第二个高斯过程G2预测噪声水平r_i，采用边缘变分方法估计异方差高斯过程回归的预测分布，称为变分异方差高斯过程；

3.5)若G3不收敛，令G1等于G3并返回3.2)，否则输出变分参数Λ和模型超参数

进一步地，在异方差高斯过程回归模型中，所述3.1)中输入向量X_*的预测分布p(y_*|X_*,X,y)服从高斯分布，即：

p(y_*|X_*,X,y)～N(m_*,∑_*)

其中，m_*和∑_*分别表示y_*的期望和协方差，N表示正态分布，样本集的大小，p(y_*|X_*,X,y)为输入向量X_*的预测分布。

进一步地，所述异方差高斯过程回归模型相应的变分，即变分异方差高斯过程回归模型为：

F(q(f),q(z))＝log p(y)-KL(q(f)q(z)||p(f,z|y))

式中，KL(.||.)表示kullback-leibler散度，q(f)和q(z)为任意的变概率密度，p(y)表示异方差高斯过程回归的边际对数似然，p(f,z|y)表示潜在变量f(x)和误差噪声z(x)的联合概率分布。

进一步地，将q(z)的分布限制为多元正态分布，即q(z)＝N(z|u,∑)，则异方差高斯过程回归模型的边缘变分MV边界可以写为：

式中，K_z表示过程z的核矩阵，其超参数为θ_z，0为N维全零矩阵，由于∫N(y|f,R)N(f|0,K_f)df＝N(y|0,K_f+R),其中R为主对角元素满足[R]_ii＝exp([μ]_i-[∑]_ii/2)的对角矩阵，则上式可进一步简化为：

式中，tr(·)表示矩阵的迹。

进一步地，所述步骤五中四项统计指标计算公式如下：

式中，

和

分别表示第i个样本的模拟值和实测值；

和

分别表示样本数据集的平均模拟值和平均实测值；N表示样本集的大小。

有益效果：

1)本发明采用完全集成的自适应噪声经验模态分解(CEEMDAN)将非线性风速时间序列分解为一组规律性更强的子序列，可以提高风速预测精度。基于CEEMDAN的预测模型可以在很大程度上提高单一机器学习模型的学习能力。

2)本发明采用变分异方差高斯过程回归模型对风速时间序列进行预测，与神经网络和支持向量机模型相比，变分异方差高斯过程回归模型的预测能力更强，且变分异方差高斯过程回归模型性能优于标准高斯过程回归模型，能够获得更高的预测精度。

附图说明

图1为本发明提供的基于变分异方差高斯过程回归的短期风速预测模型流程图；

图2为本发明实施例中样本数据集一的预测结果对比图；

图3为本发明实施例中样本数据集二的预测结果对比图；

图4为本发明实施例中样本数据集三的预测结果对比图。

具体实施方式

下面结合附图对本发明进行详细的介绍。

本发明以美国SotaventoGalicia风场在2018年每10min记录一次的风速实测数据为实施例，进行实例仿真，以验证本发明的效果。图1为本发明提供的基于变分异方差高斯过程回归的短期风速预测模型流程图，实施步骤如下：

步骤一：选取美国SotaventoGalicia风场每10min记录一次的风速实测数据作为样本数据，共选取3个数据集，每个数据集包含1008个样本数据点。根据历史实测风速数据建立风速时间序列，通过具适应性噪声的完全集合经验模态分解CEEMDAN将其分解为K个本证模态函数分量IMF_k(k＝1,2,3,...,K)和一个残差变量R_k。

步骤二：对于步骤一中每个本征模态函数分量IMF_k(k＝1,2,3,...,K)和残差R_k，计算其偏自相关函数(PACF)值，选择95％置信水平下显著的滞后时间序列作为输入变量，目标向量为输出变量，并将全部样本数据的前70％作为训练样本，剩余30％作为检验样本。

步骤三：对于步骤一中每个本征模态函数分量IMF_k(k＝1,2,3,...,K)和残差R_k，建立变分异方差高斯过程回归模型，输出变分参数Λ和模型超参数

由于变分异方差高斯过程回归是异方差高斯过程回归的扩展，因此首先介绍异方差高斯过程回归的基本理论。与高斯过程回归不同的是，异方差高斯过程G3由两个方差相同的高斯过程组成：第一个高斯过程模型称为G1，对无噪声目标输出进行建模；第二个高斯过程模型称为G2，对数噪声进行建模，即，z(x_i)＝lg(r(x_i))。该z过程由协方差函数k_z控制，其超参数为θ_z。

变分异方差高斯过程回归的步骤如下：

1)对于给定的训练集D＝{(x_i,y_i)(i＝1,2,...,n)}，其中x_i表示d维输入向量，y_i表示对应的输出向量，n表示训练样本的数量，应用第一个高斯过程G1预测输出向量y。

GPR的建模原理如下：

假设输入向量和目标输出满足以下关系式：

y_i＝f(x_i)+ε_i(1)

式中，f(x_i)表示潜在函数，其先验概率分布p(f(x_i))为高斯分布；ε_i为独立的随机变量，可表示为

2)根据G1估计训练数据的噪声水平，即z_i＝lg{var[y_i,G(x_i,X,y)]}，形成新的数据集D'＝{(x_i,z_i)(i＝1,2,...,n)}。

3)运用第二个高斯过程G2对数据集D′建模。

4)运用组合高斯过程，即异方差高斯过程G3对数据集D建模，其中，用模型G2预测噪声水平r_i，采用边缘变分(Marginalized Variational，MV)估计异方差高斯过程的预测分布，具体过程如下：

与高斯过程回归不同，异方差高斯过程回归的方差

不是常数，可以表示为

从高斯过程理论可知，给定n_*个测试数据样本D_*＝{(x_*i,y_*i)(i＝1,2,...,n_*)}，在异方差高斯过程回归模型中，输入向量X_*的预测分布p(y_*|X_*,X,y)服从高斯分布，可以表示为：

p(y_*|X_*,X,y)～N(m_*,∑_*) (2)

式中，m_*和∑_*分别表示y_*的期望和协方差。m_*和∑_*可以表示为：

m_*＝E[y_*]＝K_*(K+Ω)^-1y (3)

∑_*＝var[y_*]＝K_**+Ω_*-K_*(K+Ω)^-1K_* ^T (4)

式中，K∈R^n×n，K_ij＝k(x_i,x_j)，

K_*ij＝k(x_*i,x_j)，

K_**ij＝k(x_*i,x_*j)，Ω＝diag(r)，r＝(r(x₁),r(x₂),...,r(x_n))^T，Ω_*＝diag(r_*)，r_*＝(r(x_*1),r(x_*2),...,r(x_*n))^T；k(·)为核函数，可以表示为：

式中，k(·)为自动相关性确定平方指数(ARDSE)协方差函数函数。l_d(d＝1,2,...,D)和σ_f是模型的超参数，可表示为θ_f＝{σ_f,l₁,l₂,...,l_D}。

采用异方差高斯过程G3对数据集D建模后，输入向量X_*的预测分布p(y_*|X_*,X,y)可以转换为以下等式：

p(y_*|X_*,X,y)＝∫∫p(y_*|X_*,X,y,z,z_*)×p(z,z_*|X_*,X,y)dzdz (6)

由于式(6)右边的第二项，积分过程不能用解析的方法计算，使得p(y_*|X_*,X,y)的求解困难。因此，通常使用后验概率密度最大(MAP)求解的p(y_*|X_*,X,y)近似预测分布:

式中，

采用MAP方法能够近似估计异方差高斯过程回归的预测分布，但是MAP估计方法不能将所有的潜在变量综合起来，容易出现过拟合现象。为了克服异方差高斯过程回归的过拟合问题引入边缘变分估计异方差高斯过程回归的预测分布，称为变分异方差高斯过程回归。

对于式(2)表示的异方差高斯过程回归模型，其相应的变分近似可以定义为：

F(q(f),q(z))＝logp(y)-KL(q(f)q(z)||p(f,z|y)) (8)

由式(8)可知，对于任意的变概率密度q(f)和q(z)，异方差高斯过程回归的边际对数似然logp(y)都以F为下边界，即logp(y)≥F(q(f),q(z))。由于logp(y)的值与变分密度无关，为了获得最大的F(q(f),q(z))，需要使得KL(q(f)q(z)||p(f,z|y))的值最小，进而获得KL散度意义上的异方差高斯过程回归后验的最佳近似。由于F依赖于两个N维的变分分布q(f)和q(z)。可以通过消除F对q(f)的依赖性，进而得到一个更简单、更为紧凑的边界，称为边际变分(MarginalizedVariational，MV)边界。对于一个给定的q(z),为了消除F对q(f)的依赖性，需要计算一个新的分布q^*(f)代替q(f)，使得F(q(f),q(z))最大。q^*(f)可以通过变分贝叶斯理理论计算得到：

式中，Z(q(z))表示归一化因子。将q^*(f)代入式(8)中，可以得到异方差高斯过程回归的MV边界：

F(q(f))＝log Z(q(z))-KL(q(f)||p(z)) (10)

由于F(q(f),q(z))≤F(q^*(f),q(z))≤log p(y),式(10)表示的MV边界为异方差高斯过程回归的边界对数似然的下边界。

将q(z)的分布限制为多元正态分布，即q(z)＝N(z|u,∑)，则异方差高斯过程回归的MV边界可以写为：

式中，K_z表示过程z的核矩阵，其超参数为θ_z，0为N维全零矩阵。由于∫N(y|f,R)N(f|0,K_f)df＝N(y|0,K_f+R),其中R为主对角元素满足[R]_ii＝exp([μ]_i-[∑]_ii/2)的对角矩阵，则式(11)可以进一步简化为：

式中，tr(·)表示矩阵的迹。

由式(12)可知，异方差高斯过程回归的MV边界依赖于n+n(n+1)/2个自由变分参数，即自由参数的数量是观测样本数量的二次方。为了简化式(12)所表示的最优化问题并降低计算的复杂度，将变分参数(μ，∑)和超参数

同时优化以减少参数之间的相互作用，进一步获得只依赖于n个自由变分参数的等效边界。与GPR模型类似，对任意局部或全局最优值有，

记Λ为某个半正定对角矩阵，式(13)可以重写为：

由式(14)可知，对任意局部或全局最优值，(μ，∑)都依赖于共同的对角矩阵Λ。因此，可以根据等式(14)重新参数化(μ，∑)，此时，Λ的n个正对角元素为唯一的自由变分参数。

针对自由变分参数Λ，可以通过最小化近似和精确后验分布之间的KL发散实现下边界F(μ(Λ)，∑(Λ))＝F(Λ)的最大化。同时，针对超参数θ，也可以采用观测值边际似然的type-II最大似然(MLII)求得F的最大值。整个优化过程是非线性的且可以通过共轭梯度法完成，进而可以用解析的方式求得(Λ,θ)。由以上分析可知，变分异方差高斯过程回归的整个训练时间为O(N³)，和标准高斯过程回归的训练时间一样，因此能够快速准确地推断变分异方差高斯过程回归模型的边界。

5)若G3不收敛，令G1等于G3并返回上述步骤2)；否则输出变分参数Λ和模型超参数

步骤四：将步骤三所得的最优变分参数Λ和最优模型超参数

以及检验样本代入变分异方差高斯过程回归预测模型，即公示(8)，得到检验阶段的预测值，并对所有子模态的预测结果进行求和重构，得到原始风速时间序列的预测值。

步骤五：为了评价所提出的风速预测模型的性能，计算了预测值和实测值的均方根误差(RMSE)、平均绝对误差(MAE)、相关系数(R)和平均绝对百分比误差(MAPE)。四项统计指标计算公式如下：

式中，

和

分别表示第i个样本的模拟值和实测值；

和

采用本发明所提基于变分高斯过程回归的短期风速预测模型(CVHGPR)对3个样本风速时间序列进行预测，为了验证本发明的有效性，将其与极限学习机模型(ELM)、标准高斯过程回归模型(GPR)和不进行数据预处理的变分高斯过程回归(VHGPR)进行对比，数据集一、二、三检验期的预测结果误差指标统计值分别如表1、2、3所示：

表1 样本数据集一的预测结果误差统计

表2 样本数据集二的预测结果误差统计

表3 样本数据集三的预测结果误差统计

检验期预测结果和实测结果对比图分别如图2、3、4所示。由表可知，本发明所提预测方法预测效果优于其它三种对照模型，证明了本发明所提技术的优越性。

上述实施方式只为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡根据本发明精神实质所做的等效变换或修饰，都应涵盖在本发明的保护范围之内。

Claims

1.一种基于变分异方差高斯过程回归的短期风速预测方法，其特征在于，包括如下步骤：

步骤一：收集风电场历史实测短期风速数据，根据历史实测风速数据建立风速时间序列，通过完全集合经验模态分解将其分解为K个本征模态函数分量IMF_k(k＝1,2,3,...,K)和一个残差变量R_k；

步骤三：对于步骤一中每个本征模态函数分量IMF_k(k＝1,2,3,...,K)和残差R_k，建立变分异方差高斯过程回归模型，输出最优变分参数Λ和模型最优超参数

所述变分异方差高斯过程回归模型是引入边缘变分方法估计异方差高斯过程回归的预测分布，所述异方差高斯过程记为G3，其由两个方差相同的高斯过程组成：第一个高斯过程模型称为G1，对无噪声目标输出进行建模；第二个高斯过程模型称为G2，对数噪声进行建模，即，z(x_i)＝lg(r(x_i))，该z过程由协方差函数k_z控制，其超参数为θ_z；所述变分异方差高斯过程回归的步骤如下：