2016年第32卷第35期中国电机工程学报、

第32卷第35期中国电机工程学报、，01．32NO．35Dec。1520122012年12月15日oftheCSEE◎2012Chin．Soc．for 107ProceedingsElec．Eng．614 文章编号：0258—8013(2012)35—0107-08中图分类号：TM文献标志码：A 学科分类号：470．40采用贝叶斯一克里金一卡尔曼模型的多风电场风速短期预测卿湘运，杨富文，王行愚(华东理工大学自动化系，上海市徐汇区200237)Short·TermWindForWindFarmsSpeedForecastingMultipleUsingBayesianModeKrigedkalmanQrNGXiangyun，YANGFuwen，WANGXingyuChinaofScienceandofAutomation，East(DepartmentUniversityTechnology,XuhuiDistrict，Shanghai200237，China) ABSTRACT：Accurateshort-termwindis 计算速度更有效的变分贝叶斯方法来逼近推断和学习模型speedforecasting criticaltoreliableandsecure参数。

在公开的多风电场数据集上评估提前1h的风速预测powersystemoperations． Traditionaldidnottakeaccountofthe 性能，与持续预测算法进行比较的结果显示了该文提出的方forecastingapproachesin spatialneighborhoodsinformation．However,wind法在均方根误差评价指标上的改善。speedswindfarmsarecorrelatedbothintimeand multiplespace，This关键词：风电场；．短期风速预测；克里金一卡尔曼滤波；变awindintroducedshort-term paperspeedforecastingapproach分贝叶斯；时空模型；概率图模型modelinwhich structure usingBayesianKrigedKalman spatial 、ⅣithfunctionswasmodeledaprincipalKrigingbyBayesian0 引言 hierarchicalstructureandtimewasdynamiccomponent modeledastatemethod 将风能整合到电力系统中需要高质量的风力byspaceprocess．VariationalBayesian was toinferenceandlearnof 发电功率预测…。

风力发电不确定性的可靠估计有appliedapproximatelyparameters thewasaefficientdeterministic model，whichcomputationally利于电力系统的稳定，降低风电并网对电网的影aheadhourwasevaluatedon approach．One forecasting响，同时准确的风电功率预测能带来很大的经济收availablewinddataof windfarms．The publiclymultiple益。根据美国国家再生能源实验室和GE公司共同 resultstOweretocomparedpersistenceforecastingapproach showitsintermsofrootmean errors． 撰写的一份研究报告声【2】，美国西部电力协调委员improvementsquareelectricity KEY会(western coordinatingWORDS：windwindfarms；short—termspeedforecast；Kalman中当风电占有率为24％时，如果风电提前ld的预 Krigedfiltering；variationalBayes；spafiotemporalmodel model；probabilisticgraphical测错误能改善10％，则可以节约1亿美元的运行成 摘要：精确的短期风速预测对可靠安全的电力系统运行很重本；如果其预测错误能改善24％，则可节约1．95 要。

传统的预测方法没有考虑空间相邻风电场的信息。然而，亿美元的运行成本，因此花费较小的预测成本即可 多个风电场的风速在时间和空间上是相关的。该文给出了一取得大的经济收益。风电输出功率预测方法一般有 个采用贝叶斯一克里金一卡尔曼模型的短期风速预测方法。两类：一类是直接根据风力发电量来预测；另一类 由主克里金函数构成的空域结构使用贝叶斯层次结构进行是根据风电场风速的方向和风速通过非线性映射 建模，同时应用状态空间模型对时域动态性进行建模。采用来预测输出功率，其映射关系经常依赖风机设计。从电力系统运行来看[3】，后一类方法更可取，因为基金项目：中央高校基础研究基金项目(WHlll4059，WJ0913001)； 国家自然科学基金项目61174064)；国家重点基础研究发展相邻的风电场虽然安装的风机不一致，但是风速却 计划项目(973计划)(2012CB720502)．Funds有很强的相关性，在局部区域风速相同，各风电场FundamentalResearch fortheCentralUniversitiesNationalNaturalScience (WHlll4059，WJ0913001)；ProjectSupposedby 可根据各自的风速发电量映射关系预测输出功率。

FoundationofNationalBasicResearchChina61174064)；The 因此，本文主要研究风电场风速预测。 ProgramofChina(973Program)(2012CB720502)． 108中 国 电机工程学报第32卷电力系统运行为了平衡电力需求、供应之间的或者有些数据由于设备的原因有很明显的错误，需 矛盾和取得最小的运行代价分为不同的时间尺度要剔除掉，因此进行模型学习的数据集里经常存在 调度方法，有提前1d的经济调度(机组组合)、短期 一些缺失数据。 调度和实时自动发电控制…。相应的长期风速预测本文针对上述两个问题，提出基于贝叶斯一克 即提前至少一天的预测方法一般是基于物理方法里金一卡尔曼(Bayesian 和数值天气预报，短期风速预测则更准确，且为经电场风速短期预测方法，给出预测模型和学习方 济调度提供非常有用的信息。因此，本文主要研究法，对缺失的观测数据在学习模型参数时进行适当 短期风速预测问题，特别关注提前1h的风速预测。的处理。文献[15]率先给出克里金一卡尔曼滤波方风电场短期风速预测目前已有很多研究工作，法，在对污染和全球气候监控进行时空建模时，结 有基于神经网络、支持向量机等智能方法15-8]，也有合了两个常用的方法，即空域统计中的克里金 基于自回归、自回归滑动平均和马尔可夫过程等基 于统计模型的方法[9。

11，文献[12】则对国际上当前的模型的卡尔曼滤波方法。文献【16]则给出了贝叶斯 风力发电预测方法作了较全面的综述和评论。一些框架下的克里金一卡尔曼模型，并利用马尔可夫链chain 智能方法能取得很好的预测结果，但是对其预测机Monte蒙特卡洛(Markov 理和结果没有很好的解释，且非常依赖模型的参数学习模型参数。本文对文献【16]的贝叶斯一克里金一 整定，而基于统计模型的方法则能对预测结果给出卡尔曼模型进行了适当的修改，简化了概率模型的 合理的解释，对预测结果给出不确定性的度量方层次，增强了状态转移矩阵的灵活性，并且利用计 法。精确的点估计是风速预测的主要动机，统计模算非常有效的变分贝叶斯(variational 型不仅能利用期望值进行点估计，而且能为电力运法学习模型参数，给出了基于贝叶斯一克里金一卡尔 行和调度提供更多的关于预测结果不确定方面的曼模型的多风电场风速滚动优化和预测算法，在两 决策信息[3]o因此本文主要研究建立统计模型进行个公开的多风电场风速数据集上验证了模型的有 风速预测的方法。效性，对比持续预测法在预测精度方面有所改善。然而，目前已有的风电场风速预测研究有如下1 风电场风速的贝叶斯一克里金一卡尔曼模型 两个重要问题没有考虑或较少涉及到：I)绝大部分预测模型只考虑单个风电场风速 1．1风电场风速数据说明 预测在时间上的相关性，没有考虑多个风电场同时设有n个风电场s，，i=1…．，n，每个风电场的经 预测时在空间上的相关性。

从物理机制上看，由于 空气压力的作用接近地球表面的大气水平运动形丁个时刻的时空风速数据，设为Z(si，O，t=l，．．．，T。 成了风，覆盖了很大一片地区。在这样的空间区间1．2多个风电场风速数据的经验变异图 里不同位置的风速正相关且具有近似的特征[3】。空 1．2．1经验变异图 间距离较近的风电场风速应比较接近，空间距离远首先利用经验变异图分析空问分布的多个风 的风电场风速差异相应地较大。所以完全可以用时电场的风速数据，观察风速的空间相关性。经验 空模型对多个风电场的风速预测进行建模，既考虑变异图是一个探索性统计工具，其主要目的是显 预测在时间上的动态性，也考虑空间位置的相关示数据的空间变异。 性。空间邻近风电场的信息能提高本风电场风速预对空时风速数据Z(si，力，简单地利用一阶差分 测的精确程度。这也是本研究的主要动机和创新法去掉时域趋势分量，即：(1) 点。最新文献【13】也对短期风速预测在时空域的相w(s，，f)=Z(sf，f+1)一Z(s，，f) 关性进行了初步研究，主要探讨了风速预报误差在i=1，．．．，n；t=l，．．．，卜1 时空域的传播特性，采用的是经典的时间序列分析对一个平稳过程脚，力，传统的变异图定义为 方法；文献[14]研究了风速数据的时空统计特性。

(2) 这些研究工作主要是对风速进行分区域预测，但没以d)=研{形(Jl，t)一W(s2，f))‘】／2 有对空问相关性进行显式建模。式中d为空间位置sl和s2的距离。为简便，按文2)在收集风速数据时，有些数据不能得到， 献【16】所述，采用经验变异图公式(3)进行分析： 第35期卿湘运等：采用贝叶斯一克里金一卡尔曼模型的多风电场风速短期预测109的空间相关性。因此在建立多个风电场的风速预测触)2赤著№，f)叫％f))2(3)模型时可以考虑风电场空间位置的相关性，建立一 1．2．2测地线距离计算方法个时空数据模型将更能反映空间和时间的相关性。将空间位置s。、S：的经度和纬度值分别转为弧度 用于地统计学的克里金模型就是这样一个表征空 锄、pl和晚、易，则空间位置S1、观的测地线距离为 间相关性的模型。d=Rcos4(B) (4)而PNW风电场数据不仅存在大量的缺失数据，而且风电场在空间位置上有聚集特性，即大部 式中：R为地球半径：R=6371km；B=sin(60sin(021+分风电场聚集在一个区域内，还有几个风电场也构 cos(鼠)cos(包)cos(嗡一改)。成一个聚集，但与前面一个聚集相隔较远，因此只 1．2．3 多个风电场风速的经验变异图实例分析选取n=19个风电场的风电数据来绘制经验变异图state，以公开数据集美国威斯康星卅I(Wisconsin及指数拟合模型，如图2所示，时间点为2006年2northwest， WI)22个风电场和太平洋西北部(pacific月4日至2010年2月3日，时间间隔也为lh。

从 PNW)28个风电场的风速数据为例【17】，详细数据描述图2可以看出，太平洋西北部风电场比威斯康星州 和处理方法见实验部分。由于有大量缺失数据，WI的风电场风速空间相关性要弱。 风电场去掉5个风电场的数据，只保留n=17个风 电场数据来绘制经验变异图，时间点为2006年1吒 月1日至2009年12月31日，时间间隔为1h，即童 每隔1h记录风速数据。根据各风电场地理位置的亩咪 经度和纬度及风速数据，得到的经验变异图如图1制 所示(已将几个明显的离群点去掉)。其变异函数用 指数模型进行拟合：g(d)=Nu+s(1一exp(一3d／r))(5)图2 PNW差分数据的变异图 式中：块金值(nugget)Nu是观测误差和距离间隔很ofthedifferenceddatainPNWFig．2Variogram 小情况下变异强度的组合；S为从非零值达到较稳1．3贝叶斯一克里金一卡尔曼模型 定的常数，称之为基台值(silll，曲线开始呈水平方1．3．1 贝叶斯一克里金一卡尔曼概率层次模型 向，曲线的水平方向说明在其延伸的范围内数据对空时风电场风速数据Z(sf，吩，设在时刻t的，z 点没有空间相关性；r为变程(range)，描述了与空 间相关的差异怎样随距离变化，在变程范围内距为空间函数的随机时变线性组合： 离越近的点具有更相近的特征，因此是最重要的p(Zf)=N(Ha，，∑，)(6) 部分。

式中p(·)为随机变量的概率密度函数；脚，司为‘-、均值为∥协方差为三的正态分布概率密度函数；日攀为nxm维装载矩阵，为模型的空域分量，在本节g￥囤的结尾处阐述其构造方法。给定空间协方差函数呲暑，空间过程的无偏线性预测称之为克里金楸(Kriging)。测地线距离d／kma，为mxl矢量，表示各个风电场风速内在的时变动态性。表征空问相关的协方差矩阵暑为图1 WI差分数据的变异图inⅥ1ofthedifferenceddataFig．1Variogram五=矿／《(7)利用MATLAB的拟合函数“丘t”对WI风电场其中矩阵y的每个元素Ki： 风速的经验变异图进行非线性最小平方拟合，得到1’” Nu=O．7％=石隶忑旯如巧(弛)，202，r=182．5，拟合结果说明在空2“111(神…’扩1…～””055芦=0．1 间距离小于182．5km的范围内风速数据具有较强旯>0，r≥1(8) 110中 国电机工程学报第32卷区间参数1瓜控制当距离增加时相关性衰减的 经度和纬度值。如果考虑坐标的平方函数，则可增 速度，K“·)表示阶次为K的第二类修正贝塞尔函数。 加(经度值2，经度值×纬度值，纬度值2)3列。

F矩阵对参数丐建立一个先验概率模型：也称之为趋势场。剩下的m吲列则选择为空域方向，称为主场，p(《)～garnma(afr，％) (9) 假设空间协方差函数暑及脚均是非奇异的，计算： 式中gamma(a，61为参数为口和b的gamma分布概(14)A=(F2，F)_1F2：1 率密度函数。因此在关于Z的概率模型里，己假定空间位置B=三_1一三_1FA(15) 的相关函数对曰矩阵进行奇异值分解：B=UEU’(16)a(sf，J，)=cov(Z(sf，f)，Z(sj,f))(10)其特征矢量按特征值非降的方向排序，矩阵u和E 是平稳的，属于Matem族分布函数。都是排序后得到的矩阵，根据矩阵理论，至少存在曰动态时间趋势a。为状态空间模型：个特征值为0的特征向量，含，，z1个主克里金函数(11)p(a，)～N(Pat_1'(Q)叫)的主场则可表示为(用MATLAB的表示方法)： 应用卡尔曼滤波方法可以递归地估计状态变量值。(17)对此状态空间模型的参数也建立先验概率因此矩阵日的构造方法为H=[FpK]。 模型：2模型的变分贝叶斯学习(12)p(Q)～Wishart(2％，2b,)p(ao)～N(O，e，) (13)给定n个风电场丁个连续等间隔时刻数据zf，要根据这些观测数据学习模型参数。

传统的利用马尔 式中Wishart(a，6)为参数为a和b的Wishart先验概 率分布，P为状态转移矩阵。科夫链蒙特卡洛抽样算法计算复杂量大，需要保存很多中间抽样结果，因此提出利用确定性的变分贝叶斯因此，如式(6卜_(13)所建立的模型，组合了贝 叶斯层次模型、克里金方法及卡尔曼滤波，称之为学习方法得到模型中各随机变量的后验分布。 多风电场风速预测的贝叶斯一克里金一卡尔曼模型，变分贝叶斯学习的思想就是利用一簇简单 其概率图模型表示如图3所示，其中方框表示观测 的易处理的联合分布印@；299去逼近后验分布 变量或超参数，圆形表示随机变量，箭头表示了变p(x；互们08]，使得目标函数式(18)最大： 量之间的依赖关系。。g悸；卅式中295}日口为逼近和先验模型的参数；Z为观测数据：p@五功印@；z，嘞9(刁为观测数据与模型变量的联合似然。常见的逼近方法为均值场逼近方法，近后验分布独立。由于此模型变量的先验分布选取的均是共轭先验分布，则能很容易地写出其逼近后验分布。图3贝叶斯一克里金一卡尔曼概率图模型考虑到存在缺失值，当利用卡尔曼滤波作前向ProbabHisticforFig．3graphicalrepresentation 预测时，数据做如下处理【l州： 1．3．2空域回归矩阵日的构造方法1)将zf中所有缺失的数据设定为O；下面讲述日矩阵的构造方法。

日矩阵乘以动2)将日矩阵相应的行全部设为0。 态分量a，则表示空域回归的时变线性组合。口矩阵 因此日矩阵成为一个时变的矩阵凰。 的首先q列设为JF．，分别表示为坐标的常量、线性类似期望最大化算法，迭代计算以最大化目标 和平方函数。例如设q=3，坐标维数为2，则F的第函数(18)，得到模型中各随机变量麓的后验分布 1列元素全设为1，第2、3列则为空间坐标值，如 第35期卿湘运等：采用贝叶斯一克里金一卡尔曼模型的多风电场风速短期预测111 所示：(口，at’_1)=E(a。口21(35)Zl：丁)=ct删r+fi，扩五21lrq(r2)～gamma(tt，利用第k次迭代卡尔曼平滑计算的结果作为第r．b,)西‘=af，+砌／2斛1次迭代卡尔曼滤波计算的初始值：(19)群“=‰(36)t=％+去∑仃[y～(ZtZ：一k+l2Zt(口；)点¨+日f(口，口：)石1)]∑=q，(37)l(《)=aL／K4)在求得以上的期望值后，最大化目标函数(18)由于存在缺失数据，严格意义上来说zr要分为 可以得到矩阵P的估计： 两部分，一部分是观测到的数据，另一部分由于缺rr-1(38) 失数据的原因，需要根据其分布推出期望值，计算P=(∑(吼口2，))(∑(q彳>)一t=2t=l 较复杂，而且对结果没有太大的影响，所以采用处理论上目标函数(18)在迭代计算过程中是递增 理后的五值。

的，直至收敛。但是由于上述计算过程中涉及大量矩g(岛)～Wishart(,i日，％)阵求逆运算，影响计算速度和数值精度，因此可以采a口22a口+T-。1用简单的处理办法，即迭代一定的次数算法结束。r(20)毛=2％+∑((口f彰)一P(at一。口：)一3基于贝叶斯克里金卡尔曼模型的风电场t=2(q《1)P7+P(口，一lat7_1)P短期风速预测时间，) 风速预测算法(岛)=石。霹13．1风电场预测算法3)状态变量a，的逼近后验期望及高阶矩可用已有的风速预测研究对风速数据进行统计分 采用常规的卡尔曼滤波和平滑算法得到：析表明：未来的风速数据与相邻的短期历史数据①前向滤波：高度相关，长期历史观测数据包含预测的有用信am一1=Pfi¨卜l(21) 息。因此设计用前Ⅳ个时间点的风速数据来学习模型参数，然后再利用卡尔曼滤波算法对接下来(22)6t¨1=PGt-llt-1P’+五一天的风速进行提前1h的预测和滤波，实现滚动(23)Kt=GfIt-1I-I；(Ht6,_,耳+五)。1式优化和预测，如图4所示。由于采取了确定式(24)a啦=a州+K(五一Hta6tlt-1)的变分贝叶斯逼近推断方法，计算速度快，完全能满足预测时计算模型参数在时间上的要求。

Gtlf=(J—KtItt)ct¨1(25)学习模型参数的算法流程如图4所示。初始值为：allo=／1l(26) l学习模型参数N雾型Ⅳ+24Gllo=z1(27) —瓦……戛画磊I…丙甍4裂Ⅳ“s②后向平滑：49学习模型参数 N+48羹銎Ⅳ+72以=6tIfP7《-11If (28)(29)稚=‰+以(五Ⅲlr—a6t+lit)图4滚动优化与预测(30)Gtlr=q|r+以(G【+1Ir—Gt+1]f)∥horizonandFig．4Recedingoptimizationforecasting 则后验期望及高阶矩可由卡尔曼平滑结果得到：输入：给定n个风电场Ⅳ个时间间隔．1h的风(a，)=E(a短期风速预测时间，IZl：r)=五灯(31)e最大迭代次数。(32) 一；(吩口；)=f(a，a’lzl：r)=GtIr+％％初始化：设定模型参数，