基于贝叶斯方法的风力发电短期概率预测

(整期优先)网络出版时间:2020-09-28
/ 4

基于贝叶斯方法的风力发电短期概率预测

孙一睿

中国电力工程顾问集团西北电力设计院有限公司

陕西省西安市 710000

摘要

目前,风力发电在可分布式发电系统中占重要地位,给社会生产和人民生活带来了极大的经济和环境效益。然而,风能的间歇性和不稳定性使得风力发电系统难以发挥其最佳操作性能并满足市场的需求。因此,如何准确预测风力发电的发电量已成为近年来广泛关注的问题。已发表的文献多采用确定性预测的方式对风电产量进行预测,但近年来概率预测的方法已引起了广泛关注。

本文提出了一种先进的概率预测方法,用于风电产量的短期预测,使用两种韦布尔分布的混合概率函数来模拟风速的不确定性。然后,使用自回归、集成、移动平均模型的贝叶斯推断方法来确定混合韦布尔分布的参数。

关键词:风能;风力发电; 预测方法; 概率方法

  1. 引言

本文提出了一种短期概率预测的新方法,该方法可以直接对风力发电的输出进行预测。该方法利用风电有功功率与风速之间的关系,通过贝叶斯推断(BI)方法预测其概率密度函数(PDF)。

众所周知[1],基于贝叶斯方法的短期概率风速预测的两个关键步骤是:

(1)建立影响风速的不确定性因素的PDF分析表达式;

(2)定义最佳时间序列模型以排除非贝叶斯方法的先验随机参数的PDF参数。

参考文献[2,3,4,5]中介绍了一些关于以上问题的研究。在文献[3]中,作者使用高斯PDF对风速进行建模,并采用了仅包含风速的六阶自回归模型。在文献[4]中,作者采用了正态分布的混合模型来拟合失速速度的接近值,而使用韦布尔分布来拟合其剩余值。在本文中,将采用更为复杂的PDF来模拟风速的不确定性,并采用一种有效的方法来确定最佳时间序列模型,同时采用两个韦布尔分布的混合分布作为PDF分析表达式。

本文考虑采用两个韦布尔分布的混合分布,是因为在文献[6]中得出了这样的结论:仅包含两个参数的经典韦布尔分布不能模拟自然界中出现的所有风态,例如具有双峰分布的风态。因此,必须为每种风态匹配更合适的PDF,以最大程度地减少产能估算中的误差。实验证明,两种韦布尔分布的混合函数更适合用来模拟单峰和双峰风态。

本文的研究目的是:

(1)提出一种新型基于贝叶斯公式的风电短期预测方法;

(2)在贝叶斯方法的框架内采用新的概率函数和改进的时间序列模型;

(3)对新贝叶斯方法、传统贝叶斯方法和参考预测方法(概率持续性方法)的性能进行比较,总结该方法的优缺点。

  1. 基于贝叶斯方法的风力发电量预测

本文研究的是贝叶斯方法在预测风力发电系统产生的有功功率的PDF(概率密度函数)方面的应用,着重研究风力发电有功功率与风速之间的关系,并在蒙特卡洛模拟方法的框架内应用所选定的关系。

2.1风能有功功率与风速之间关系

一般情况下,风电有功功率不仅取决于风速,还取决于气象变量,如风向、温度、局部空气密度和降水量。在许多情况下,还需要对整个风电场的风能进行预测,由于风力发电场中发电机可能具有不同的切入速度和额定速度,这将导致确定性功率曲线的选择变得极为复杂。发电设备的增加或者对已有设备进行维护都有可能引起风电场的变化 [7-9]。

然而,在相关文献中经常采用确定性功率曲线[7-9]。实际上,厂家提供的有关功率曲线通常是基于空气密度恒定、环境变量取标准值的假设。

考虑到风电的有功功率是一个随机变量,它不仅仅取决于风速还取决于其他显性变量(例如风向或空气密度),这将导致需要估算的参数数量的增加,从而增加贝叶斯推断的复杂性。

时间范围h内,有功功率5f71aae655ace_html_dd13af52ca1c63e8.gif和风速5f71aae655ace_html_e5b2291993846698.gif在之间的以下解析关系可以写成:

5f71aae655ace_html_1f20b3a2ea9512ce.gif (1)

其中:

5f71aae655ace_html_3abf8dee7e8a0cf4.gif是由线性函数、抛物线函数或三次函数近似的非线性函数;

Wci,Wr和Wco分别是风力发电机组的切入、额定和截止特性值。

2.2. 风速PDF解析表达式的选择

众所周知,风速通常使用韦尔布分布(WB)进行建模,如文献[18]中所述:

5f71aae655ace_html_f5529c3627999932.gif (2)

其中:

5f71aae655ace_html_31efc7b3f20cc82f.gif是比例参数;

5f71aae655ace_html_bb6359566ae2ac8c.gif是形状参数。

比例参数5f71aae655ace_html_31efc7b3f20cc82f.gif可以根据风速分布的平均值5f71aae655ace_html_f9f18dd797b28e92.gif表示,其关系式如下:

5f71aae655ace_html_7e7435a8f867a153.gif (3)

因此,可以将公式(2)中的PDF视为平均值和形状因子的函数。

文献[6]得出的结论是,两个参数的韦布尔分布具有一系列的优点,但韦尔布PDF不能代表自然界中遇到的所有风况,比如双峰分布的风况,而两种韦布尔分布的混合将更好地模拟此类情况。

双组分混合的韦布尔分布(MWB)取决于五个参数(5f71aae655ace_html_e5b2291993846698.gif5f71aae655ace_html_60a00bd35abd0a38.gif5f71aae655ace_html_ca4aa4ebbf77bdcd.gif5f71aae655ace_html_b399341fc0935eb1.gif5f71aae655ace_html_ac9b54b7632e98f6.gif),并由下式给出:

5f71aae655ace_html_23dcb6a0d6899db.gif

5f71aae655ace_html_19a2445caa86cfe2.gif(4)

其中:5f71aae655ace_html_196df943aa0dee6f.gif

公式(4)中的比例参数5f71aae655ace_html_60a00bd35abd0a38.gif可以根据以下关系用风速分部的平均值5f71aae655ace_html_9e7bfba04ee3deca.gif和其他参数5f71aae655ace_html_748088350f91af59.gif表示[6]:

5f71aae655ace_html_69e472db0dd0b5dc.gif (5)

由等式(4)和(5)的分析结果可知,在时间范围h之内通过求解风速平均值

5f71aae655ace_html_9e7bfba04ee3deca.gif和估计参数5f71aae655ace_html_748088350f91af59.gif可以准确地预测概率密度函数5f71aae655ace_html_47521674df173d3e.gif。 在本文中,假设参数5f71aae655ace_html_748088350f91af59.gif是贝叶斯方法的先验随机参数,而5f71aae655ace_html_9e7bfba04ee3deca.gif的值将使用下一小节提到的自回归移动平均值(ARMA)和ARIMA时间序列模型来估算。

2.3.回归移动平均值(ARMA)和ARIMA时间序列模型

随机变量5f71aae655ace_html_4fb54ac5dc1c301d.gif的通用ARMA表达式可以表示为:

5f71aae655ace_html_1a1d1c649d0bd21b.gif (6)

其中:

B是后向移位算子,由5f71aae655ace_html_e3427b851385eedb.gif确定;

5f71aae655ace_html_b8c9bc6f11759719.gif是由5f71aae655ace_html_3e3812983772d236.gif定义的p阶平稳自回归算子;

5f71aae655ace_html_20a2fdecfef5558e.gif是一个常数项;

5f71aae655ace_html_693db1b27f96d80c.gif是q阶的移动平均算子; 5f71aae655ace_html_d977efad9e504b4.gif

5f71aae655ace_html_54eb7139c5dfc412.gif是时间t处的白噪声,具有零均值和恒定方差5f71aae655ace_html_d01234c31ccf7251.gif

根据5f71aae655ace_html_4fb54ac5dc1c301d.gif5f71aae655ace_html_54eb7139c5dfc412.gif的过去值扩展方程式(6),我们获得以下形式的差分方程式:

5f71aae655ace_html_7055c4dc4ae36ad2.gif (7)

因此,通过确定阶次(p,q)和p + q + 2未知参数5f71aae655ace_html_20a2fdecfef5558e.gif5f71aae655ace_html_3ae2f2872814ea09.gif,…,5f71aae655ace_html_53daba3bd98b45c2.gif,5f71aae655ace_html_54a130af99b1b518.gif5f71aae655ace_html_2dc8c0e6fea360d4.gif,5f71aae655ace_html_4e0835d4bd8f81bd.gif可以确定ARMA模型。

ARMA模型以数学方式表示线性、平稳的随机过程,但是当用来模拟非平稳过程时,此模型的性能通常较差。为了更准确地用数学方式表示这些时间序列,必须使用ARMA的扩展版本。

此类模型属于ARIMA类,对于通用随机变量5f71aae655ace_html_4fb54ac5dc1c301d.gif,可以表示为:

5f71aae655ace_html_3b4c151f0b6744cf.gif (8)

其中:

5f71aae655ace_html_e06ba1520ec3ba30.gif是由5f71aae655ace_html_961baeafaa0b95b9.gif定义的后向差分算子。应注意,多项式Φ(B)必须满足上述平稳条件。

根据5f71aae655ace_html_4fb54ac5dc1c301d.gif5f71aae655ace_html_54eb7139c5dfc412.gif的过去值扩展方程式(8),我们得出以下形式的差分方程式:

5f71aae655ace_html_82e532c97445e7fb.gif(9)

其中:系数5f71aae655ace_html_40288e1624752497.gif,5f71aae655ace_html_81e6c6a89cd14cb.gif是算子5f71aae655ace_html_aa0c81aac8ac90ef.gif的系数。实际上,多项式5f71aae655ace_html_b8c9bc6f11759719.gif可以看成是两部分的共同作用,即具有d=1解的多项式5f71aae655ace_html_ef2cfaffaf2b0eb9.gif和表示上述平稳性要求的多项式5f71aae655ace_html_f482e6de6f5b47eb.gif

因此,通过确定模型阶数(p,d,q)和求解未知参5f71aae655ace_html_5684b8df387e8847.gif,…,5f71aae655ace_html_81e6c6a89cd14cb.gif,5f71aae655ace_html_54a130af99b1b518.gif5f71aae655ace_html_2dc8c0e6fea360d4.gif,5f71aae655ace_html_7aff5633737ee8d8.gif可以确定ARIMA模型。在本文中,我们在样本自相关函数5f71aae655ace_html_7395ce6093c96780.gif的基础上使用Box-Jenkins方法[19],该方法是对不同滞后5f71aae655ace_html_8b6d4776f1ec1eb4.gif下的理论自相关函数5f71aae655ace_html_fd06a55c795a5906.gif的估计:

5f71aae655ace_html_ba71c39344d2ea5a.gif (10)

其中5f71aae655ace_html_13b0bb1b48fb615.gif5f71aae655ace_html_334d1cad3ef98490.gif分别是随机变量5f71aae655ace_html_4fb54ac5dc1c301d.gif的理论均值和理论方差。由于时间序列始终由有限数量的样本N组成,因此只能提供5f71aae655ace_html_fd06a55c795a5906.gif的估计值5f71aae655ace_html_8c5328f97930a84e.gif,如下所示:

5f71aae655ace_html_8df5c9593d88b212.gif (11)

其中5f71aae655ace_html_b9d327cd634ba473.gif是时间序列的样本均值。

Box-Jenkins方法的第一步是利用自相关函数确定差异度d。 实际上,对于平稳的时间序列,样本自相关函数5f71aae655ace_html_8c5328f97930a84e.gif对于中等滞后量5f71aae655ace_html_8b6d4776f1ec1eb4.gif迅速衰减为零;而对于非平稳的时间序列,自相关函数5f71aae655ace_html_8c5328f97930a84e.gif的下降非常缓慢,即使存在较大的滞后也不趋于零。这表明:

如果样本自相关函数5f71aae655ace_html_8c5328f97930a84e.gif随着5f71aae655ace_html_8b6d4776f1ec1eb4.gif的值增加而迅速减小,则该时间序列可以用平稳模型表示;如果样本的自相关函数5f71aae655ace_html_8c5328f97930a84e.gif随着5f71aae655ace_html_8b6d4776f1ec1eb4.gif值增加没有迅速减小,则此随机过程在5f71aae655ace_html_4fb54ac5dc1c301d.gif处是不稳定的,但5f71aae655ace_html_cd826b7ae6876926.gif的情况下在5f71aae655ace_html_e26ded06206440ef.gif处是稳定的。

一旦确定了微分阶数d,则经过适当微分的时间序列5f71aae655ace_html_8560367a7ba18e80.gif将显示平稳过程的特征。因此,可以通过(p,q)阶ARMA模型对其进行建模。以这种方式建立了时间序列5f71aae655ace_html_2234edfe0a8b2f9e.gif后,代表5f71aae655ace_html_2234edfe0a8b2f9e.gif的ARMA(p,q)模型和代表原始时间序列5f71aae655ace_html_4fb54ac5dc1c301d.gif的ARIMA(p,d,q)模型具有相同的阶数p,q。因此,在Box-Jenkins方法的第二步中,可以通过研究差分时间序列5f71aae655ace_html_2234edfe0a8b2f9e.gif的阶数来确定原始ARIMA模型的阶数。

研究结果表明,微分时间序列5f71aae655ace_html_2234edfe0a8b2f9e.gif的自相关函数5f71aae655ace_html_8c5328f97930a84e.gif随着阶数(p,q)的取值不同有不同的表现。表1列出了最常见时间序列取值:

表1 5f71aae655ace_html_eca48f5bd26c6b0e.gif与阶数(p,q)的关系

ARIMA模型的阶数

(1,d,0)

(0,d,1)

(2,d,0)

(0,d,2)

(1,d,1)

5f71aae655ace_html_8c5b9e500cee2230.gif呈指数下降

5f71aae655ace_html_6854e4822a09bfca.gif是唯一可知的非零项

5f71aae655ace_html_8c5b9e500cee2230.gif是指数函数或正弦波的混合

5f71aae655ace_html_50107cacf9f7e2d6.gif是唯一可知的非零项

5f71aae655ace_html_8c5b9e500cee2230.gif5f71aae655ace_html_6854e4822a09bfca.gif之后呈指数下降

一旦确定了ARIMA模型的三个阶数(p,d,q),就可以使用估算的方式来确定等式(9)中p+q+2个未知参数5f71aae655ace_html_5684b8df387e8847.gif,…,5f71aae655ace_html_53daba3bd98b45c2.gif,5f71aae655ace_html_54a130af99b1b518.gif5f71aae655ace_html_2dc8c0e6fea360d4.gif,5f71aae655ace_html_7aff5633737ee8d8.gif的值,进而可以确定时间序列模型。

2.4参数5f71aae655ace_html_55d43154a70dca8b.gif5f71aae655ace_html_ac9b54b7632e98f6.gif的概率密度函数估算

一旦如第2.3节所述确定了时间范围h处的分布平均值,则必须求得方程(4)中的其余参数5f71aae655ace_html_748088350f91af59.gif。贝叶斯推断先根据每个参数的已知(或假设)先验概率分布推断一组观察值,然后对这些参数进行概率估计,最后确定它们的联合后验概率分布。

初始集合5f71aae655ace_html_d7cfca462d4dd4a1.gif是由初始时间h-k处观测到的M个风速测量值组成,然后选择参数的先验分布。设5f71aae655ace_html_1b228b9fb85a3b9c.gif为通用参数,其在时间范围h内的先验分布必须首先确定,先验分布的参数通常称为超参数。在本文中,将初始估算值5f71aae655ace_html_1b228b9fb85a3b9c.gif应用估算算法代入风速5f71aae655ace_html_12ade6d6ea178360.gif集合的观测值对各参数5f71aae655ace_html_f6cc522719d96e58.gif在进行估计。然后,假定每个估算的结果等于相应的先验信息高斯分布的平均值,即 5f71aae655ace_html_f2f527d8c41d3b02.gif

一旦设置了先验分布,BI便可以通过贝叶斯定理的扩展,在给定测量值5f71aae655ace_html_12ade6d6ea178360.gif的情况下估算联合后验分布5f71aae655ace_html_72857d036755a1a9.gif。然而,无法通过解析地提供5f71aae655ace_html_72857d036755a1a9.gif的闭式表达式,但是与5f71aae655ace_html_72857d036755a1a9.gif成正比的非标准化后验分布5f71aae655ace_html_3d0570330a41ba44.gif的表达式对于参数的概率估计而言已足够。

我们可以通过以下等式计算随机参数的非标准化后验分布

5f71aae655ace_html_3d0570330a41ba44.gif

5f71aae655ace_html_f5c1c4af4a498020.gif (12)

其中:

5f71aae655ace_html_bdfab8af495920ad.gif是似然函数;

5f71aae655ace_html_df7ca2ee1156e2c.gif是向量5f71aae655ace_html_1b228b9fb85a3b9c.gif的第j个先验随机参数的先验分布。公式(12)中的似然函数5f71aae655ace_html_bdfab8af495920ad.gif由下式给出:

5f71aae655ace_html_aa4e934ee43a6d82.gif (13)

其中,5f71aae655ace_html_fd2571179329058e.gif对应于由公式(5)求出的参数5f71aae655ace_html_1b228b9fb85a3b9c.gif5f71aae655ace_html_2d82f00eaac817fa.gif,其中5f71aae655ace_html_2d82f00eaac817fa.gif是在给定包含在集合5f71aae655ace_html_1c62c243ffe6a956.gif中的原始风速5f71aae655ace_html_b9b3f71d612cd22a.gif的p+d值的前提下,从选定的ARIMA(p,d,q)模型得出的最小均方误差预测值。

似然函数的显式表达式如下:

5f71aae655ace_html_5a715ff0c66f054.gif (14)

可以通过相关文献中使用的各种贝叶斯方法求得方程式(12)的数值。 本文采用基于Metropolis-Hasting算法的蒙特卡洛马尔可夫链模拟方法,通过对未归一化后验分布5f71aae655ace_html_108d1073ff5f971e.gif的求解,可以得出5f71aae655ace_html_1b228b9fb85a3b9c.gif中参数后验分布的样本。

2.5概率密度函数5f71aae655ace_html_2232307ac961c42c.gif的样本估计

2.4节中求得的时间范围h内参数5f71aae655ace_html_f3139b0098f4f380.gif的后验分布样本和第2.3节中求得的平均值5f71aae655ace_html_b3d6ccba44807c17.gif可以一起使用,通过公式(5)获取参数5f71aae655ace_html_60a00bd35abd0a38.gif的样本。然后,可以从估算分布5f71aae655ace_html_d57eaf4d449e77e6.gif中获取风速样本5f71aae655ace_html_805f13b4fc4bd89d.gif。然后可以在蒙特卡洛程序中使用5f71aae655ace_html_805f13b4fc4bd89d.gif的样本,以从等式(1)获得5f71aae655ace_html_27f4d4e5e0a864f2.gif的样本,从而进行风有功功率的概率估计。

  1. 结论

在本文中,我们提出了一种新型概率方法对风力发电进行预测。该方法基于贝叶斯理论,适用于在风速随时间变化很大的情况下来预测风能。该过程是通过使用复杂的概率分布(混合韦布尔分布)表征风速来实现的,此种方法尤其适合用来表示单峰和双峰风态。

将该方法的求解结果与两种已广泛使用的概率预测方法求得的结果进行比较发现该方法对风电的短期概率预测十分有效,在点值和分布预测方面都比参考方法更加准确,且此方法在预报的可靠性方面也有显著的提升。

我们得出的结论是,本文所述贝叶斯方法最适合于表示单峰和双峰风况。在风速长期以双峰分布为特征的特殊情况下,预测工作还会有进一步的改善和提高。

未来的研究工作将集中在基于样本自相关函数的Box-Jenkins方法的应用上,以预测光伏发电量。另外,由于本文使用了确定性功率曲线,下一步将根据风向、温度、局部空气密度和降水等气象变量对功率曲线的概率进行研究。

参考文献

[1] Gelman, A.; Carlin, J.B.; Stern, H.S.; Rubin, D.B. Bayesian Data Analysis; Chapman & Hall:London, UK, 1995.

[2] Bracale, A.; Caramia, P.; Carpinelli, G.; Di Fazio, A.R.; Ferruzzi, G. A Bayesian method for short-term probabilistic forecasting of photovoltaic generation in smart grid operation and control. Energies 2013, 6, 733–747.

[3] Miranda, M.S.; Dunn, R.W. One-Hour-Ahead Wind Speed Prediction using a Bayesian methodology.In Proceedings of the IEEE Power Engineering Society General Meeting, Montreal, QC, Canada,18–22 June 2006.

[4] Zhang, J.; Pu, J.; McCalley, J.D.; Stern, H.; Gallus, W.A., Jr. A Bayesian approach for short-term transmission line thermal overload risk assessment. IEEE Trans. Power Deliv. 2002, 17, 770–778.

[5] Bracale, A.; Caramia, P.; Carpinelli, G.; Di Fazio, A.R.; Varilone, P. A Bayesian-based approachfor a short-term steady-state forecast of a smart grid. IEEE Trans. Smart Grid 2013, 4, 1760–1771.

[6] Carta, J.A.; Ramırez, P.; Velazquez, S. A review of wind speed probability distributions used in wind energy analysis. Case studies in the Canary islands. Renew. Sustain. Energy Rev. 2009, 13, 933–955.

[7] Taylor, J.W.; McSharry, P.E.; Buizza, R. Wind power density forecasting using ensemble predictions and time series models.

IEEE Trans. Energy Convers. 2009, 24, 775–782.

[8] Hering, A.S.; Genton, M. Powering up with space-time wind forecasting. J. Am. Stat. Assoc. 2009,105, 92–104.

[9] Gneiting, T. Quantiles as optimal point forecasts. Int. J. Forecast. 2011, 27, 197–207.