基于时空协方差函数的风能场景生成方法与应用

彭星皓, 李艳婷

(上海交通大学 机械与动力工程学院, 上海 200240)

为实现碳达峰、碳中和的目标,满足日益增长的能源需求,可再生能源已被广泛应用于能源系统和市场.风力发电(简称风电)是一种在世界范围内储备丰富的可再生能源,已成为许多地区越来越重要的发电来源.截至2021年底,全球风电累计装机容量达到840 GW,并且仍处于快速上升阶段[1].对风力或风电功率时间序列的点和概率性预测是风电预测系统的重要组成部分.然而,对于经济调度这类随机优化问题,点和概率性预测并不能提供充分参考.近年来,场景生成法已成为风电预测的一个新兴研究领域.场景生成指通过构建数学模型,对未来一段时间内风电场的风速或功率做出预测,生成多个符合风电场时间序列和时空特征的场景.生成场景作为经济调度、机组组合等优化决策问题的输入,对合理地求解调度方案、进行电力系统运营相关的决策具有重要意义.在风能场景预测中,一个场景通常指风电场在一个时间段内的风能功率时间序列.

总体而言,场景生成方法包含预测、抽样以及优化方法等.预测方法中常用的预测模型如自回归和移动平均(Auto Regressive and Moving Average, ARMA)[2-3]是用于刻画序列线性相关性的参数化工具,在预测单一风电场时往往具有较好效果,但难以刻画多个地点间的非线性关系;另一预测方法是机器学习,其包含的神经网络模型如生成对抗网络(Generative Adversarial Network, GAN)[4-6]和Wasserstein 生成对抗网络(Wasserstein Generative Adversarial Network, WGAN)[7]等方法无需考虑原始数据分布的统计假设,理论上数据样本越多,生成场景的质量越高,但对于训练样本的规模和质量依赖较大.抽样方法中Copula建模[8-10]、拉丁超立方抽样[11]、蒙特卡洛马尔可夫链抽样[12]较为常见,其中Copula是一种普遍使用的概率分布建模方法,可对单一风电场的时间相关性、对多个风电场的空间相关性进行建模,也可构建风电和其他可再生能源之间的关系[13],但受构建的概率密度是否可逆、维度灾难等问题影响,此方法计算负担大.优化抽样方法是场景生成的延伸,它通过矩匹配[14]或距离匹配[15]来削减场景,应用于场景生成后的进一步优化.

现阶段,大规模风力发电源于广阔区域内多个风电场,需要对多个风电场总功率场景做出预测,因此把握风电场之间的空间相关性和风电场内的时间相关性非常重要;同时,具体的场景生成依赖于功率预测值和真实值之间的联合分布,因此准确构建预测值和真实值的条件联合分布是另一重点.当前主流场景生成方法对这两个重点有不同处理方式.Wang等[15]利用核密度估计(Kernel Density Estimation, KDE)法建立功率真实值和预测值间的条件概率函数,并使用R-vine Copula构建风电场间的空间相关性,最后通过条件概率逆变换来得到生成场景.但在R-vine Copula构建过程中,需要依次估计不同风电场间条件概率分布的Copula函数,受风电场数量影响,计算过程繁琐复杂,且只考虑了风电场间的空间相关性,忽视了风电场内时间相关性对建模的影响.Ma等[16]将风能预测值按大小分层,根据不同区间预测功率对应的实际功率构建条件概率分布,并使用指数型协方差函数表示风能时间相关性,但指数协方差函数只考虑了时间相关性,不能反映风电场的空间相关性,只适用于单一或少数风电场.Tan等[17]使用高斯混合模型估计功率真实值和预测值间的条件联合概率分布函数,但对数据分布的要求是混合高斯分布,对预测值和真实值的联合分布的估计不够准确.Deng等[18]利用C-vine Copula构建风电场内部的空间相关性结构,但空间相关性只体现于两个风电场之间,且并未生成具体的每日风电出力场景.

针对目前场景生成方法中存在只考虑风电功率的时间相关性或只考虑空间相关性以及条件概率函数估计不够准确等问题,本文创新在于利用结合经验估计方法的时空协方差函数,准确刻画了多个风电场功率序列的时间和空间相关性特征,同时使用Pair Copula模型准确拟合功率预测值和实际值之间的联合分布,生成具体每日风电出力场景,并通过指标进行评估.此外,为验证生成的风能功率场景在电力系统中的运用价值,介绍一种两阶段的电力系统经济调度模型,并考察模型可行性与可靠性,进一步验证方法生成场景质量.

1 数据集介绍

使用美国中东部地区2012—2013年中25个风电场的Wind Toolkit数据集[19],风电场标号以及空间分布如图1所示,其中坐标点上方数字为风电场标号.该数据集包含了每个风电场的风速和功率数据,每个风电场额定功率为10～16 MW.功率数据集由100 m轮毂高度处的风力数据和适合站点位置的风力发电机功率曲线创建,考虑风力发电机间的唤醒效应以估算每个站点产生的功率,分辨率为5 min.为减少计算量,取分辨率为1 h.利用长短期记忆(Long Short-Term Memory, LSTM)神经网络的方法进行训练和预测得到预测数据集,均方根误差为1.71 MW.

2 基于时空协方差函数场景生成方法

总体方法框架如图2所示,首先介绍场景生成方法所需数学知识,然后具体介绍时空协方差函数构建方法即Pair Copula构建方法,最后给出场景生成具体步骤.

图2 本文方法总体框架图Fig.2 Overall framework of method in this paper

2.1 预备知识

2.1.1时空协方差函数 设x(h1,u1),x(h2,u2)为发生在时刻u1,u2和位置h1,h2的两个时空过程,其相关性可由时空协方差函数刻画.Gneiting[20]提出不可分割的时空协方差函数反映时空过程的联系.时空协方差函数一般如下式所示:

(1)

式中:u,h分别为时间间隔和空间距离;σ2为整体数据的方差;a,c为时间和空间距离的尺度参数;α,γ分别为时间协方差和空间协方差的平滑参数;β∈[0, 1]为时空协方差的关联,当β=0时,表示时间协方差和空间协方差分割;参数τ≥β表示对时间协方差做出调整.不可分割的时空协方差函数优势在于可以针对使用的数据集进行灵活的参数调整,从而提高描述时空过程相关性的准确性.

2.1.2Kendall秩相关 Kendall秩[21]涉及和谐对概念,设w,z为两个长度为n的观测序列,对于分别来自两个序列的两对观测值(wi,zi)和(wj,zj),若(wi-zi)(wj-zj)≥0,则称这是一个和谐对,反之为非和谐对,Kendall秩相关系数的计算公式如下:

(2)

式中:n1,n2分别为和谐对和非和谐对数量.Kendall秩相关系数反映了变量间的协变关系,在没有线性假设和数据非正态分布的情况下适用于刻画风电场功率时空相关性情景.而最常见的皮尔逊相关系数则通常应用于刻画线性相关的变量相关性.

2.1.3Copula理论 在Copula理论中,Sklar定理[22]解释了多元分布函数、Copula函数和边际分布函数之间的关系.假设一个d维随机变量x=[x1x2…xd]T的边际累积分布函数(Cumulative Distribution Function, CDF)为F1,F2, …,Fd,其联合累积分布函数为F,那么若所有的边际函数都连续,则F可以被一个特定的Copula函数C定义:

F(x)=C(F1(x1),F2(x2), …,Fd(xd))

(3)

概率密度函数(Probability Density Function,PDF)表示为

f(x)=

(4)

式中:c(·)为Copula的概率密度函数;fi(·)为边际概率密度函数.

2.2 时空协方差函数的参数估计

在计算时空协方差矩阵之前,需要对数据的时空平稳性进行检验,以保证时空协方差函数效果.首先通过通过增广Dickey-Fuller(Augmented Dickey-Fuller, ADF)检验验证时间序列平稳性.对于空间平稳性,在缺少空间平稳性假设时需要进行多维标度变换[23-24]使序列空间平稳.根据文献[23],本文用新坐标取代原始经纬度坐标表征风电场间空间相关性,在变换后的平面中,空间相关性仅通过距离表示,其空间相关性与原始平面中一致,且变换坐标后的数据具有时空平稳性假设,为后续时空协方差函数估计做好铺垫.坐标变换前后风电场间距离和相关系数的散点图如图3所示,反映变换前空间距离(h)和空间相关性的负相关关系并不明显,但经过多维标度变换后空间距离和空间相关性呈现明显负相关关系,便于时空协方差函数更好地捕捉时空相关性.

图3 坐标变换前后风电场间距离和相关性系数的散点图Fig.3 Scatter plots of spatial distances between wind farms and correlation coefficients at different coordinates

图4 不同方法的等高线图Fig.4 Contour plots of different methods

根据文献[20]参数设置,在协方差函数参数中,σ2,γ,τ均取为1.在进行ξ={a,c,α,β}估计时,通过对比经验的时空相关性等高线图调整参数.在调整后的参数附近,对目标函数:

(5)

进行优化,使用网格搜索法对参数进行微调.最终求得时空协方差函数等高线图如图4(b)所示.

将不同时间间隔和空间距离代入协方差函数,可以得到表示25个风电场24 h实际功率的时空相关性协方差矩阵Σ,矩阵维度为600×600,热力图如图5所示,其中位于(24j1+t1, 24j2+t2)的小格代表j1风电场第t1小时的功率与j2风电场第t2小时的功率的协方差(0≤t1,t2≤23 h, 1≤j1,j2≤25).当j1=j2时,热力图表示同一风电场的协方差矩阵,可以看出在同一风电场中,时间越接近则协方差越大,时间相关性越强;当j1≠j2时,热力图显示不同风电场间的协方差,因为地理位置上25个风电场主要分成两个簇(见图1),所以协方差矩阵的热力图具有两个深色区块为前16个风电场和后9个风电场,在同一区块中的风电场距离接近,协方差较大,时空相关性较强.

2.3 基于Pair Copula模型联合概率分布估计

时空协方差函数捕捉时空过程相关性,但不能直接生成具体场景.风电功率具体场景的生成依赖于实际功率和预测功率之间的关系.条件的联合概率分布反映了已知功率预测值时实际功率取值的分布情况,因此快速准确估计功率预测值和实际值间的联合概率分布对后续的场景生成尤为重要.

根据Copula理论,利用Copula函数构建实际功率和预测功率联合概率分布.由于只有实际功率和预测功率两组数据,所以使用连接两个变量的Pair Copula函数来连接两者的分布最合适.常见的Pair Copula函数有Gaussian Copula、t-Copula以及Archimedean Copula族中的Frank Copula、Gumbel Copula、Clayton Copula等,该类Copula能够对包括非对称依赖性在内的复杂依赖性结构进行建模.

以68736号风电场2012年的实际功率数据和预测功率数据为例,其联合分布形式如图6所示.可以看到功率的实际值和预测值在极值附近联合分布密集,其他区域分布相对稀疏,且大致对称,符合Frank Copula概率密度函数的分布特点,因此采用Frank Copula构建联合概率分布.

图6 68736号风电场2012年风电实际功率和预测功率的散点图Fig.6 Scatter plot of actual and predicted wind power of No. 68736 wind farm in 2012

设v1,v2为两个随机变量,θ为待求参数,Frank Copula的函数表达式如下:

(6)

在具体计算过程中,需要先将数据转化为累计分布函数值.设pr和pf分别为某一风电场j每小时的功率实际值和预测值,v1和v2为对应的经验累积分布函数(Empirical Cumulative Distribution Function, ECDF)值,则v1=FECD(pr),v2=FECD(pf),联合概率密度函数fj=cj(v1,v2),累积概率分布函数的边际概率为均匀分布,概率密度为1.Frank Copula只含一个参数,在使用极大似然估计时较为方便.图7为θ=23.055时,68736号风电场数据估计的Frank Copula概率密度函数.

图7 68736号风电场数据的Frank Copula概率密度函数Fig.7 Frank Copula probability density function plot of the data of No. 68736 wind farm

为验证Frank Copula的准确性,选取同为单参数的Gumbel Copula和Clayton Copula进行对比,指标为相对于经验Copula的平方距离(d2),如表1所示.可知Frank Copula与经验Copula距离最小,准确性最高.同时由该Copula函数随机抽样生成真实值与预测值联合分布模拟值,得到对应分位数-分位数图如图8所示,可知模拟值与原始值的分布情况基本一致,证明了Frank Copula的准确性.

表1 不同的Copula与经验Copula的平方距离Tab.1 Squared distance of different Copulas from empirical Copulas

图8 真实值和预测值与其模拟值的分位数-分位数图Fig.8 Quantile-quantile plots of measured values and predicted values with their simulated data

2.4 场景生成步骤

对于典型多维高斯过程,可以使用随机生成多维高斯分布随机数模拟.通常,风电功率的分布是非正态的、风电功率的相关性也是非线性的[25],但可以通过基于时空协方差函数的高斯过程生成风电功率的条件概率值,最后通过逆变换求得具体场景值,具体步骤如下:

(7)

(3) 对于j=1, 2, …,Nj,t=1, 2, …, 24 h,n=1, 2, …,Ns,重复步骤(2),得到每个风电场24 h的所有场景.将一个风电场每小时的对应场景值组合,得到每日的一个场景,不同风电场的场景值累加得到累积场景X.

3 场景评估方法与算例结果对比分析

首先介绍几种常用的场景评估指标,然后利用数据集提供的历史数据生成7个风电场90 d的场景,最后评估与分析不同场景生成方法所生成的场景.

3.1 场景评估方法介绍

(8)

3.1.2变异函数得分 变异函数得分(Variogram Scoring, VS)既反映生成场景相对于真实值的精度,又可以用于展示不同场景生成方法对时空相关性的捕捉能力[27],变异函数得分越小表示该方法捕捉时空相关性能力越强.p阶VS计算真实值与生成场景值的变差函数之间的差异定义为

(9)

本文使用的阶数p设为0.5,权重wij设为1.

3.1.3覆盖概率 覆盖概率(Coverage Probability, CP)表示生成的场景值包含真实值的概率,反映了场景值能否捕捉极端情况的发生的能力,表达式如下:

(10)

3.1.4场景区间宽度 场景区间宽度(bSI)反映场景值相对于真实值的集中程度,可以用于评估场景生成方法捕捉不同风电场间的空间相关性的能力.场景区间宽度由平均区间宽度(bAI)和带系数λ的平均区间偏差(bAW)组成,bAI表示了场景值高估或低估真实值的误差,在经济调度问题中,越低的bAI意味着功率预留量越低,为应对电力系统运行中因为预测误差而分配的向上和向下功率储备所造成的浪费,bAI越小越好,相关表达式为

(11)

3.2 场景生成结果与对比分析

图9 不同场景生成数量下的风电功率波动情况的核密度估计Fig.9 Kernel density estimation of wind power fluctuations with different scenario numbers

同时采取其他3个场景生成方法作为基准进行对比,本文提出的方法记作M1,文献[15]提出的Vine-Copula以及加权核密度估计拟合联合分布记作M2,文献[16]提出的根据预测功率分别构建条件概率分布记作M3,文献[17]提出的协方差函数和高斯混合模型估计联合概率分布记作M4,但未根据功率波动情况进行聚类与场景削减.使用Python语言编程,计算在配置为Intel(R) Core(TM) 2.60 GHz CPU和16 GB内存的个人计算机上进行.M1～M4生成未来24 h每小时100个场景的平均用时分别为74、65、34、52 s.由于是对未来1 d的场景进行预测,1 min左右的场景生成时间相对于实际24 h调度周期在可接受范围,4种场景生成方法均具备时间效率,所以不在计算时间上比较优劣.

表2 4种方法生成的场景在各评价指标下的对比Tab.2 Comparison of scenes generated by four methods at each evaluation index

为直观展示本文方法的效果,随机选取某日(2013年1月21日)生成的场景进行绘图.如图10所示,本方法生成的场景能够准确捕捉到实际功率(P′)的变动趋势,距离真实值的分布相较于其他方法更紧密.在生成的每日100个场景的基础上,本文通过场景削减方法提取出10个典型的场景,其中聚类数量可以选择.具体方法为通过K-means聚类算法,将100个场景聚集成10个类别,每个类别取类内均值作为一个典型场景,每个类别的概率等于该类别所包含的场景数量的1/100.2013年1月21日的10个典型场景如图11所示,其中P为每个典型场景的概率,Ptot为风电场总功率.

图10 2013年1月21日4种方法生成的风电场总功率场景Fig.10 Total wind power scenarios generated by four methods on January 21, 2013

图11 2013年1月21日10个典型的风电场总功率场景及其概率Fig.11 10 typical total wind power scenarios and their probabilities on January 21, 2013

4 场景生成的应用

介绍生成的风能场景在现实的工业背景中的应用方法,包括电力系统经济调度中的一种两阶段优化模型,构建相关的机组组合混合整数规划模型,求解以本方法生成的场景在某日调度中的优化结果并进行分析.以4种方法生成的场景分别作为输入,对90 d调度计划的平均成本进行计算与对比.

4.1 两阶段联合调度问题介绍

电力系统经济调度指在满足安全和电能质量的前提下,合理利用能源和设备,以最低的发电成本或燃料费用保证对用户可靠供电的一种调度方法.传统的电力系统依赖燃煤或燃气等发电单位来平衡供给和需求.然而,当风电使用水平提高时,由于风电的错峰性和不可调节性,所以系统需要更多备用容量或灵活资源调度来应对不确定性.从发电方面来看,灵活发电技术[28]、区域储备优化[29-30]和储能[31]是常见选择.本文主要考虑灵活发电技术中对可调控机组的调度问题,在该问题中,系统运营商协调调度常规燃煤或燃气机组以应对风电波动.通常以整体运营成本最小为目标,以可调控的燃煤或燃气机组等发电机组每小时的出力大小、机组的启停状态、风电削减量以及缺负荷量等可以人为调控的参数为决策变量,以机组运行、供需平衡等限制为约束进行优化求解.主要的决策变量为发电机组的出力情况,且决策变量既包含发电单位的启停状态这类整数变量又包含出力大小这类连续型变量,因此此类问题被称作电力系统的机组组合问题,需通过建立混合整数规划模型求解.

一种常见的联合调度的总体框架如图12所示.图中调度考虑了用户响应,即通过阶梯电价策略调控用户用电需求,产生虚拟发电效果.出于问题简化的考虑忽视了用户响应,考虑调度中心对传统热电站的实时出力和预留功率调度,以及不同风电场景下的风电削减和非自愿用户用电削减问题,这种联合调度可被分为两阶段优化:第一阶段为日前调度,调度员根据风电功率和负载的预测值确定备用容量和发电机组启停状态以及大致的出力,其中备用容量指未来1 d实际场景的波动在每一时段预留的发电功率容量;第二阶段对应于未来1 d日内操作期间的部署,即根据实际用电量和实际风电场发电量实时调整热电站的出力.在日前调度中,调度员根据当前的预测值做出第一阶段决策,因此将其定义为此时此地(Here and Now, HN)阶段.第二阶段当风电实际出力、实际负荷等不确定的参数给出时,调度员在原有的调度计划基础上做出必要调整,因此称其为等待和观望(Wait and See, WS)阶段.

图12 电力系统两阶段联合调度框架图Fig.12 Framework of two-stage joint dispatching of power system

在一个典型的电力系统联合调度问题中,优化目标为最小化调度过程的两个阶段的总成本:

minCtotal=min(CHN+CWS)

其中:CHN为第一阶段的成本;CWS为第二阶段的成本,且

CHN=Cfuel+Csetup+CRUD

(12)

(13)

式中:Cfuel为能耗成本;Csetup为启动成本;CRUD为预留容量成本;ρn为每个场景出现的概率;CAFn为每个场景下的能耗调整成本;CLSn为非自愿功率削减成本也称需求未满足成本;CWCn为风力削减成本.需要满足的约束包括总功率平衡约束、相邻时段的功率变化范围约束、发电机组开启和关停时间约束、输电线容量限制、预留功率容量限制等约束,具体可参考文献[32].

4.2 两阶段联合调度问题求解与分析

本系统将IEEE 6总线系统的一个热电机组替换为风电场,如图13所示.其中G1、G2、G3代表3个可调控的热电机组,WP代表风电场,Load 1、Load 2、Load 3代表负载.3个热电机组的最大有功功率分别设置为250、250、200 MW,负载每小时的总功率在300～600 MW区间波动,在每个场景中,负载功率仅发生轻微波动.因为生成的场景能够合理反映风电的实际波动情况,所以将生成的大量场景作为第二日的实际风电情况输入至模型,以验证模型的可行性和可靠性.首先对单日优化调度进行计算与分析,然后对比4种场景生成方法生成的场景输入90 d调度中的求解结果.使用Gurobi 9.1.1对以上混合整数规划问题进行求解,精度设为0.1%.

表3 不同场景生成方法的优化求解结果对比Tab.3 Optimization solution comparison of different scenario generation methods 美元

图13 调整后的IEEE 6总线系统示意图Fig.13 Schematic diagram of adjusted IEEE 6-bus system

风电的一大特性是错峰性,风电功率的波峰通常出现在夜间即负荷较少的时刻,选用2013年1月21日这一典型错峰出现的日期生成场景.取每个场景出现的概率ρn=1/Ns,带入相关参数进行求解.模型总求解时间为3.42 s,证明了模型的潜在可行性.第一阶段的求解结果如图14(a)所示,第二阶段的求解结果随机抽取一个场景作为展示如图14(b)所示.由图14(a)可知,在第一阶段中需求较小时,仅需要开启一个热电机组配合风电进行供给;在需求较大时,需开启多台机组.从4号场景的WS阶段调度方案可以看出,此场景在7～11 h内,风能比预期较少,已开启的两个机组已经满功率运行,因此出现需求未满足的情况;15～20 h内此场景下风能较多,出现了风力削减的情况.

图14 电力系统两阶段调度求解结果Fig.14 Two-stage scheduling solution results of power system

4.3 不同场景生成方法的应用与对比

参照4.2节的设置,分别将4种场景生成方法在2013年1月1日至3月31日90 d内生成的场景作为两阶段调度计划的输入进行求解,并计算各项成本的均值,结果如表3所示.

可知,燃料成本和机组启停成本基本一致,原因为其仅和两阶段联合调度的第一阶段有关;但在燃料调整成本、风力削减成本、需求未满足成本等项目上,本文方法具有明显优势,原因为本文预测场景值与真实值更接近,场景的分布更紧凑,第一阶段的调度策略则更精准有效,减少了第二阶段的调整与浪费以及由预留容量造成的损失.就总成本而言,本方法相较于方法M2～M4分别减少了3.61%、3.16%和5.40%,体现了解决大规模的风电并网的电力系统调度问题的经济性.

5 结论

针对风电功率场景生成现存的时空相关性刻画不全面、联合条件概率分布描述不准确等问题,提出一种基于时空协方差函数、Pair Copula以及高斯过程概率逆变换的场景生成方法,通过在Wind Toolkit数据集上的应用以及在不同指标上的结果对比,得出以下结论:

(1) 时空协方差函数能够同时对多个风电场功率的时空相关性进行有效建模,具有运用方便灵活的特点.

(2) 在Pair Copula中,单参数的Frank Copula能够较准确地建立实际功率和预测功率之间的条件概率函数.

(3) 本方法生成的场景准确把握了实际场景的变化趋势,相较于其他方法,在多项场景评估指标上均有较好表现,表明本方法能准确抓住多个风电场功率的时空相关性,具有捕捉极端场景的能力,并能通过聚类缩减为少数几个典型场景,在风能场景预测问题上具有应用价值.

风能场景的生成最终需要应用于电力系统中,作为经济调度、机组组合等决策优化问题的输入,生成场景有助于电力系统操作员为未来的调度方案做出决策.建立了未来24 h的两阶段电力调度的混合整数规划模型.通过对比实验,验证了方法的成本优势,并证明了本文场景生成方法在风电并网的实际问题中具有潜在应用价值.出于计算方便考虑,该模型对部分情况进行简化,未来将考虑更符合实际的电力调度模型,考虑多种新能源的场景生成方法,并在具有更复杂机制的经济调度模型中应用.