1.风电功率预测物理模型
风电功率预测最先采用的是物理方法,早期的物理方法类似欧洲风图集的推理,经过近几十年的发展,物理方法获得了较大的发展。风电功率物理预测方法技术路线如图2.3所示,主要包括数值天气预报数据引入、考虑风电场区域粗糙度和地形变化的风机轮毂高度处风速、风向获取以及风速-功率转化三个技术环节。
数值天气预报数据引入环节包括:选取与风机轮毂高度临近的数值天气预报层高的风速和风向;引入数值天气预报中的气温和气压对风机的标准功率曲线进行修正。
风机轮毂高度处参考风速、风向获取环节较为复杂。由图2.3可知主要包括对数风廓线模型、地转拖曳定律、粗糙度变化模型和地形变化模型。
(1)对数风廓线模型
对数风廓线用于描述中性近地面层中风速随高度变化的情况,可求取不同高度的风速,是广泛采用的一种风廓线。中性层结下的风廓线为
图2.3 风电功率物理预测方法技术路线
式中, 为获得的指定高度 z 下的平均风速; u * 为具有速度量纲的非负常数,称为摩擦速度,近地面层摩擦速度不随高度变化; κ 为卡门常数,其数值大约在0.30~0.42之间,一般取0.40; z 0 为地表粗糙度,对应平均风速为0的高度。已知层高 z 1 处的平均风速 ,层高 z 2 处的平均风速 , z 0 可通过下式计算获得:
(2)地转拖曳定律
根据准地转近似,可将大尺度系统下由简化运动方程计算得到的地转风作为自由大气中实际风的近似,若忽略热成风的作用,可认为地转风在大尺度系统中保持不变,且维持水平匀速直线运动。地转风在大尺度系统中保持不变的特点,使得地转风可作为联系大气边界层中不同位置风速、风向的桥梁。根据地转拖曳定律,可建立地转风与近地面层特征量摩擦速度 u * 与地表粗糙度的关系公式:
式中, G 为地转风; f 为地转参数,取值与经纬度有关,中高纬度地区可设为10 -4 s -1 ;经验常数 A 和 B 依赖于大气层结稳定度,在中性层结下有 A =1.8, B =4.5。地转风与地表风的风向存在夹角,其正弦值可计算为:sin α =-( Bu * )/( κ | G |)。
若已知某未知处的风速、风向,并获得了该区域的有效粗糙度,则可根据式(2.6)求得摩擦速度 u * ,进而由式(2.8)求得地转风。在地转风大尺度范围内保持不变的假设前提下,可进一步由地转风求得风电场内任意位置的风速和风向。
(3)粗糙度变化模型
在实际中,利用地转风推导风机轮毂高度风速和风向需考虑风电场内地表粗糙度变化的影响,该过程可通过对数风廓线体现,使研究位置处的风速表现为不同对数风廓线在不同层高下的拼接,如图2.4所示。
图2.4 粗糙度变化下对数风廓线示意图
假设上游未受干扰来流经过两次粗糙度变化扰动后到达风机所在位置, x 1 和 x 2 分别为引起第1次和第2次粗糙度变化的物体到风机的距离,此时,风机轮廓线应由3个部分拼接而成,分别为对应粗糙度 z 01 、摩擦速度 u *1 的 u 1 ( z ),对应粗糙度 z 02 、摩擦速度 u *2 的 u 2 ( z ),以及对应粗糙度 z 03 、摩擦速度 u *3 的 u 3 ( z )。根据实验观测与仿真分析,流经变化粗糙度的下风向风廓线公式可描述为
式中, z 02 为风机位置的粗糙度; z 01 为初始来流的粗糙度; ,其中 u *1 和 u *2 分别为 z 01 和 z 02 对应的摩擦速度; h 为内边界层高度,由下式计算确定:
式中, z ' 0 =max( z 01 , z 02 ), x 为粗糙度变化位置与研究位置的距离。粗糙度与摩擦速度的关系为
根据以上结果,进而根据风电机组地理位置、内边界层高度 h 以及风电机组轮毂高度,可由式(2.8)得到粗糙度变化影响下风电机组轮毂高度处风速。
(4)地形变化模型
地形变化主要影响边界层风速变化,其对边界层的影响可认为是相对于上风向未受扰对数风廓线的一个小扰动,并可将地形扰动下的受扰边界层分为内外两层。
对于边界层外层,其在地形扰动下的流场变化可按照势流理论求解;对于边界层内层,其变化为非线性惯性力项与湍流摩擦力共同作用的结果,因此流场扰动随高度按对数风廓线变化。
风速-功率转化环节通过风机功率曲线实现,风机功率曲线有标准功率曲线和实际功率曲线两种。如果风电场具备经专业机构测试的实际功率曲线,则应采用实际功率曲线,如不存在实际功率曲线,则需采用理论功率曲线。理论功率曲线是在标准空气密度下给出的,需根据预测的气温和气压对其进行修正,修正公式如下式所示:
式中, ρ 为空气密度; M 为气体摩尔质量; P 为气压; R 为比例常数; T 为气体的热力学温度。
2.光伏功率预测物理模型
光伏功率预测的物理方法即根据光伏发电的物理原理和光电转换效率来确定光伏电站的输出功率。首先利用数值天气预报等提供的辐照度预报值,经过预处理后结合光伏电站的地理位置及光伏电池板倾角等信息,采用太阳位置模型得到光伏电池板接收的有效辐照强度;然后利用构建的光电转换效率模型,将光伏电池板接收的有效辐照强度转化为输出功率;最终经过全场累加获得整个光伏电站的输出功率预测结果。
光电转换效率模型是整个物理预测过程的关键,需要建立经验公式并合理地确定经验参数。光伏电池组件输出直流功率为
P ( t )= ηAI (2.13)
式中, P ( t )为输出功率; η 为光电转换效率; A 为光伏电池组件的总面积; I 为太阳有效辐照度。
光电转换效率的常用模型有以下四种,在实际预测过程中,可根据情况合理选择:
1)常系数效率模型:
η = η STC (2.14)
式中, η STC 为常量,根据材料取不同的数值。
2)考虑温度影响的模型:
η = η STC [1- β ( T -25℃)] (2.15)
式中, T 为环境温度; β 为常量,取值范围为0.003~0.005。
3)考虑太阳辐射和温度的模型:
η ( I , T )=( a 1 + a 2 I + a 3 ln( I ))[1- β ( T -25℃)] (2.16)
式中, a 1 ~ a 3 为经验参数,可由最小二乘法拟合求解。
4)考虑太阳辐射、温度、大气质量影响的模型:
式中,AM为大气质量; p 、 q 、 m 、 r 、 s 、 n 为经验参数,可通过不同工况测试得到。
1.持续法
新能源发电功率序列在短时间内往往表现出一定的时序惯性,即临近时刻的功率在数值上较为相近。基于此特性,持续法将前一时刻的量测功率直接作为未来时刻的功率预测值。持续法模型简单直接,在气象变化平稳时具有优异的预测表现,然而当气象波动较大时预测精度将受到较大的影响。
2.ARIMA模型
差分自回归移动平均(Autoregressive Integrated Moving Average,ARIMA)模型,是对自回归移动平均(Auto Regression Moving Average,ARMA)模型的改进。ARMA模型适用于平稳非白噪声序列的数据,然而在实际应用中所面临的数据大多为非平稳数据,因此需要在ARMA模型的基础上结合差分运算进行改进,即为ARIMA模型。
ARMA模型可以近似表示为
Y t =β 0 + β 1 Y t -1 + β 2 Y t -2 +…+ β p Y t - p + ε t + α 1 ε t -1 + α 2 ε t -2 +…+ α q ε t - q (2.18)
式中, β 1 、 β 2 、…、 β p , α 1 、 α 2 、…、 α q 为常数;{ ε t }为白噪声序列。称时间序列{ Y t }为服从( p , q )阶的移动平均模型,简称为ARMA( p , q ),或者记为 φ ( B ) Y t = θ ( B ) ε t ,其中 φ ( B )=1- β 1 B - β 2 B 2 -…- β p B p 称为滞后算子多项式,同理, θ ( B )=1+ a 1 B + a 2 B 2 +…+ a q B q 。
通常来说,产生时间序列的某一随机过程可以由三部分组成,如下式所示:
Y t = f t + P t + X t (2.19)
式中, f t 为趋势项; P t 为周期项; X t 为随机项。当 f t 和 P t 为变量时,表明时间序列为非平稳随机过程。
实际中存在的随机过程大多为非平稳随机过程,ARMA模型的平稳条件是方程 φ ( B )=0的根必须位于单位圆外。对于非平稳时间序列,需先利用差分方法将其平稳化。非平稳时间序列{ Y t }进行一阶有序差分后如下所示:
∇ Y t =(1- B ) Y t = Y t - Y t -1 (2.20)
经过 d 阶差分后公式为
∇ d Y t =(1- B ) d Y t (2.21)
则原ARMA模型 φ ( B ) Y t = θ ( B ) ε t 经过 d 阶差分后可表示为
φ ( B )∇ d Y t = θ ( B ) ε t (2.22)
即ARIMA模型( p , d , q )。
建立ARIMA模型的步骤包括平稳性检验、非白噪声检验、模型定阶以及最后的预测,其中模型定阶的常用方法有:残差-方差图定阶法、自相关函数和偏自相关函数定阶法、F检验定阶法以及最佳准则函数法。
1.支持向量机
支持向量机(Support Vector Machine,SVM)是一种监督学习模型,广泛应用于各种分类问题和回归问题。支持向量机的决策目标是寻找一个能最大化训练样本分类间距的超平面。当样本线性可分时,支持向量机可直接在原空间寻找最优分类超平面;当样本线性不可分时,支持向量机利用非线性映射将低维度空间的输入样本映射到高维度空间,使其变为线性可分,并引入松弛变量,进而在该高维空间中寻找最优分类超平面。相比于逻辑回归和神经网络,支持向量机在学习复杂非线性方程时提供了一种更为清晰、更加强大的思路。
给定训练样本集{( x i , y i ), i =1,2,…, N , x i ∈ R n , y i ∈ R },支持向量机的工作原理是寻找一个非线性映射 ϕ ( x ),将数据 x 映射到高维特征空间 F 中,并在该高维特征空间 F 中利用以下估计函数 f ( x )进行线性回归:
f ( x )=[ ωϕ ( x )]+ b , ϕ : R m → F , ω ∈ F (2.23)
式中, ω 为权重向量; b 为偏置值。其函数逼近问题等价于如下函数最小化:
式中, R reg [ f ]为期望风险; R emp [ f ]为经验风险,且 ; λ 为常数。
通过构造损失函数并基于结构风险最小化的思想,支持向量机构建如下优化问题来确定模型参数:
式中, C 为用来平衡模型复杂项和训练误差项的惩罚参数; ξ * i 、 ξ i 为松弛因子; ε 为不敏感损失函数。该问题可转化为以下对偶问题:
式中, a i 与 a * i 为拉格朗日乘子; K ( x i , x j )为核函数,求解式(2.26)可得到支持向量机回归函数为
根据支持向量机回归函数的性质,只有少数 a * i - a i 不为零,这些参数对应的向量成为支持向量,向量机回归函数 f ( x )完全由其决定。核函数可以用原空间中的函数实现,无需了解该非线性变换的具体形式,常用的核函数有如下几种:
1)多项式核函数:
K ( x , y )=( xy +1) d , d ≥1 (2.28)
2)径向基核函数:
3)Sigmoid核函数:
K ( x , y )=tanh[ b ( xy )+ θ ] (2.30)
2.极限学习机
极限学习机(Extreme Learning Machine,ELM)是一种单隐含层前馈神经网络模型。相比于传统的神经网络模型,极限学习机结构简单,不用设置学习率,具有学习速度快、训练误差小、泛化能力强的优点。
给定训练集样本{( x i , y i ), i =1,2,…, N , x i ∈R n , y i ∈R},ELM的隐含层单元个数为 k ,激活函数为 g (·),则ELM的输出模型为
式中, β j 为连接第 j 个隐含层节点和输出节点的权重; a j 为连接第 j 个隐含层节点和输入节点的权重向量; d j 为第 j 个隐含层节点的偏置矩阵; g (·)为激活函数。
ELM模型的待训练参数包括 a j 、 β j 、 d j ,在训练过程中满足下式:
该过程可由矩阵表示为
Hβ = Y (2.33)
从而,连接权重 β 可通过最小化2-范数 得到,即
β = H + Y (2.35)
式中, H + 为隐含层输出矩阵 H 的Moore-Penrose广义逆矩阵。
综上所述,极限学习机的具体建模步骤为
1)确定激活函数 g (·)以及隐含层单元个数 k 。
2)随机生成权重矩阵 a j 以及隐含层偏置矩阵 d j 。
3)根据已知量求得隐含层输出矩阵 H 及其摩尔-彭罗斯(Moore-Penrose)广义逆矩阵。
4)根据式(2.35)求出连接权重 β 。
3.BP神经网络
BP神经网络是一种由具有自适应性的简单单元构成的广泛并行互联的网络,它的组织结构能够模拟生物神经系统对真实世界所做出的交互反应。神经网络的一般结构包括输入层、隐藏层和输出层,各层之间都采用全连接的形式。全连接的本质为矩阵、向量相乘,仍然是线性变换的过程,为了提高神经网络的非线性拟合能力,往往在输出层中增加激活函数。以单隐含层神经网络为例,其架构可用下式表示:
H = f ( U * I + b I ) (2.36)
Y = g ( V * H + b H ) (2.37)
O = σ ( Y ) (2.38)
式中, I 为神经网络的输入向量, H 为隐藏向量; Y 为激活前输出向量; O 为输出向量; f 、 g 和 σ 为激活函数; U 和 V 为权重矩阵, b I 和 b H 为偏置向量,*为矩阵乘法算子。
BP神经网络中常用的激活函数有Sigmoid、Tanh和Relu函数:
1)Sigmoid函数:
2)Tanh函数:
3)Relu函数:
f ( x )=max(0, x )(2.41)
给定包含 N 个样本的训练集,可根据神经网络输出的预测结果 O 和真实值 Z 构建损失函数如下:
BP神经网络的训练过程是基于训练样本不断优化权重矩阵 U 和 V ,使得损失函数极小化。在BP神经网络中,式(2.36)~式(2.38)被称为神经网络的正向传播过程。通过正向传播可以得到BP神经网络的输出结果,进而得到训练的损失函数,将损失函数的原函数对权重矩阵求偏导,用来更新权重矩阵,该过程叫作神经网络的反向传播过程。反向传播的公式如下:
δ O =-( Z k - O k )× σ′ ( Y k ) (2.43)
δ H =( δ O * V )× σ′ ( Y k ) (2.44)
δ I =( δ H * U )× σ′ ( Y k ) (2.45)
式中, δ O 、 δ H 和 δ I 为损失函数对各层的梯度参数; 和 为损失函数对各层的偏导数值; σ′ ( Y k )为激活函数对激活值的偏导。进而可通过下列公式对各层权重矩阵进行更新:
式中, α 为学习率。当隐含层为多层时,仍可根据链式求导法,即按照式(2.43)~式(2.48)的思路求各层偏导,进而对权重矩阵进行更新。