CN103530531A - 一种基于极大似然估计的风电功率持续特性描述方法 - Google Patents

Abstract

本发明公开了一种基于极大似然估计的风电功率持续特性描述方法，利用该方法分析，发现风电功率的不同状态的持续时间特性与逆高斯分布相似，同时，风电场有可能长时间保持相同的出力状态，且在出力较低或较高的状态保持不变的平均时间要长于保持中等出力状态的平均时间；风电功率持续时间分布特性与逆高斯函数分布曲线大致吻合的这一发现，将有助于对风电功率在某一出力水平下的持续时间做出估计，对于风电功率的短期预测提供重要参考；本发明有助于对风电功率在某一出力水平下的持续时间做出估计，对于风电功率的短期预测提供重要参考；对含风电电力系统的规划、运行方式的安排等，都有一定的指导意义。

Description

一种基于极大似然估计的风电功率持续特性描述方法

技术领域

本发明涉及风力发电技术领域，尤其涉及一种基于极大似然估计的风电功率持续特性描述方法。

背景技术

据统计，我国新增风电装机占全球新增装机的比例从2006年的不足10%上升到2012年的35%。截至2012年6月，全国并网风电容量已经有5258万千瓦，居世界首位。由于风资源的不确定性和风电机组本身的运行特性，使风电场的输出功率具有间歇性和波动性，大规模风电的接入势必给电力系统的安全稳定运行带来困难，故加强对风电的随机特性的研究，对提高风电场出力的预测精度，进而提高电网运行水平，从而降低非可再生能源的消耗，保障电力系统安全稳定，提高电力系统经济性，减少温室气体排放具有重大意义。目前，对风电功率随机特性中的波动特性研究较多，而对于持续特性尚缺乏合理的定义与定量的分析。因此，如何定量地描述风电功率持续特性的概率分布是一个难点问题。

发明内容

本发明的目的就是为了解决上述问题，基于风电功率状态的定义，提出状态持续时间指标，利用极大似然估计法，发现风电功率状态的持续特性服从逆高斯分布，完善风电功率特性分析体系。

为了实现上述目的，本发明采用如下技术方案：

一种基于极大似然估计的风电功率持续特性描述方法，包括以下步骤：

（1）定义风电功率的状态：

将风电功率的可能取值范围离散化为若干个功率区间，每个功率区间即为风电功率的一个状态；依照定义，实测风电功率序列中的每一个数值均对应一个风电功率状态；

（2）进行风电功率状态持续时间的概率统计：

风电功率状态持续时间的统计包括两个方面：保持某个状态的时间长度和持续某一个时间长度的次数；例如：当风电功率从任意状态m，进入到状态n后，其中m≠n，开始记录风电功率保持在状态n内的时间；若风电功率经历时间T后跳出状态n，则记录状态n持续时间T一次；按照这种方法统计实测风电功率序列，可以得到风电功率在状态n下不同持续时间各自出现的次数，分别除以对应状态n下持续时间的总次数即得到该状态n下持续时间的概率，统计该状态n下所有持续时间T的概率值，则可以得到风电功率在此状态下持续时间的概率分布；

（3）选择适当的分布函数：

选择Birnbaum-Saunders分布、指数分布、逆高斯分布以及对数正态分布作为参考分布函数；

（4）利用最小二乘估计法和极大似然估计分别对选择的分布函数进行拟合；

（5）采用残差平方和指标计算所选择函数的拟合精度。

所述步骤1中，风电场的额定装机容量记为P_E，拟划分的状态总数为N，则第n个状态代表的功率区间范围设定为：

$\begin{array}{c}(P_{\mathrm{lower}}^n,P_{\mathrm{upper}}^n]\\ P_{\mathrm{lower}}^n=(n-1)\×\frac{P_E}{N},P_{\mathrm{upper}}^n=n\×\frac{P_E}{N},n=\mathrm{1,2},...,N\end{array}---\left(1\right)$

其中，Pⁿ _upper和Pⁿ _lower分别代表第n个状态所代表的风电功率范围的上下限，此外，风电场由于无风或弃风导致输出功率为零的情况较多，因此将零值专门定义为风电功率的一个状态，即0状态。

所述步骤4中，各个函数拟合方法如下：

A)Birnbaum-Saunders函数

其中，为标准正态密度函数，α、β为形状参数和尺度参数，均为未知参数。

这里采用最小二乘估计计算其未知参数:

由(2)式得

$\begin{array}{c}\frac{1}{\mathrm{\α}}(\frac{\sqrt{x}}{\sqrt{\mathrm{\β}}}-\frac{\sqrt{\mathrm{\β}}}{\sqrt{x}})~N\left(\mathrm{0,1}\right)\\ \frac{\sqrt{x}}{\sqrt{\mathrm{\β}}}-\frac{\sqrt{\mathrm{\β}}}{\sqrt{x}}~N(0,{\mathrm{\α}}^2)\end{array}---\left(3\right)$

记： $\mathrm{\ϵ}=\frac{\sqrt{x}}{\sqrt{\mathrm{\β}}}-\frac{\sqrt{\mathrm{\β}}}{\sqrt{x}}---\left(4\right)$

则： $\frac{\sqrt{x}}{\sqrt{\mathrm{\β}}}=\frac{\sqrt{\mathrm{\β}}}{\sqrt{x}}+\mathrm{\ϵ}$

其中，ε～N(0，α²)        (5)

可见，可看作一个回归模型。设来自x的样本为：x₁、x₂、x₃Lx_n，先求β的最小二乘估计。

令 $\begin{array}{c}Q\left(\mathrm{\β}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(\frac{\sqrt{x_i}}{\sqrt{\mathrm{\β}}}-\frac{\sqrt{\mathrm{\β}}}{\sqrt{x_i}})}^2\\ \frac{\mathrm{\δQ}}{\mathrm{\δ\β}}=2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(\frac{\sqrt{x_i}}{\sqrt{\mathrm{\β}}}-\frac{\sqrt{\mathrm{\β}}}{\sqrt{x_i}})(-\frac{\sqrt{x_i}}{2\mathrm{\β}\sqrt{\mathrm{\β}}}-\frac{1}{2\sqrt{\mathrm{\β}}\sqrt{x_i}})=0\end{array}---\left(6\right)$

解得 $\mathrm{\β}=\frac{1}{n}\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{x_i}}---\left(7\right)$

下面求α的点估计，记误差平方和为s_E,则由回归分析的结果，可取α²的估计为

${\mathrm{\α}}^2=\frac{S_E}{n-1}=\frac{1}{n-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(\frac{\sqrt{x_i}}{\sqrt{\mathrm{\β}}}-\frac{\sqrt{\mathrm{\β}}}{\sqrt{x_i}})}^2=\frac{2n}{n-1}(\frac{1}{n}\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{x_i}}-1)---\left(8\right)$

B)指数分布函数

$f(x;\mathrm{\λ})=\left\{\begin{array}{cc}{\mathrm{\λe}}^{-\mathrm{\λx}}& x>0\\ 0& x\≤0\end{array}\right.---\left(9\right)$

式中，λ为率参数，即每单位时间内发生某事件的概率，

这里采用极大似然估计法计算其未知参数λ：

$L\left(\mathrm{\λ}\right)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}{\mathrm{\λe}}^{-\mathrm{\λ}{x}_i}={\mathrm{\λ}}^n{e}^{-\mathrm{\λn}\stackrel{\‾}{x}}---\left(10\right)$

其中：

是样本均值

似然函数对数的导函数为：

$\frac{d}{\mathrm{d\λ}}\mathrm{InL}\left(\mathrm{\λ}\right)=\frac{d}{\mathrm{d\λ}}(\mathrm{nIn}\left(\mathrm{\λ}\right)-\mathrm{\λn}\stackrel{\‾}{x})=\frac{n}{\mathrm{\λ}}-n\stackrel{\‾}{x}---\left(11\right)$

得λ估计为：

C)逆高斯分布函数

$f(x;\mathrm{\μ},\mathrm{\λ})={\left[\frac{\mathrm{\λ}}{2\mathrm{\π}{x}^3}\right]}^{1/2}\exp \frac{-\mathrm{\λ}{(x-\mathrm{\μ})}^{\text{2}}}{2{\mathrm{\μ}}^{2}x}---\left(12\right)$

式中，μ>0为均值，λ>0为形状参数。在λ相同的情况下，μ值越大，则分布的尖峰越低。当λ趋近于无穷时，逆高斯分布逐渐趋近于正态分布；

这里采用极大似然估计法计算其未知参数μ和λ：

极大似然函数为：

$L(\mathrm{\μ},\mathrm{\λ})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}f(x_i;\mathrm{\μ},\mathrm{\λ})=\frac{{\mathrm{\μ}}^n}{{\left(2\mathrm{\π\λ}\right)}^{\frac{n}{2}}\underset{i=1}{\overset{n}{\mathrm{\Π}}}{\left(x_i\right)}^{\frac{3}{2}}}{e}^{\{-\frac{1}{2\mathrm{\λ}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\mathrm{\μ})}^2}{x_i}\}}---\left(13\right)$

两边取自然对数得：

$\mathrm{InL}(\mathrm{\μ},\mathrm{\λ})=\mathrm{nIn\μ}-\frac{n}{2}\mathrm{In}\left(2\mathrm{\π\λ}\right)-\frac{3}{2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{Inx}}_i-\frac{1}{2\mathrm{\λ}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\mathrm{\μ})}^2}{x_i}---\left(14\right)$

分别求关于μ,λ的偏导数，得似然方程组：

$\begin{array}{c}\frac{\mathrm{\δInL}(\mathrm{\μ},\mathrm{\λ})}{\mathrm{\δ\μ}}=\frac{n}{\mathrm{\μ}}+\frac{1}{\mathrm{\λ}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{x_i-\mathrm{\μ}}{x_i}=0\\ \frac{\mathrm{\δInL}(\mathrm{\μ},\mathrm{\λ})}{\mathrm{\δ\μ}}=-\frac{n}{2\mathrm{\β}}+\frac{1}{2{\mathrm{\β}}^2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\mathrm{\μ})}^2}{x_i}=0\end{array}---\left(15\right)$

解得： $\mathrm{\μ}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i,\mathrm{\β}={\left(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i\right)}^2\left(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{x_i}\right)-\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i---\left(16\right)$

D)对数正态分布

$f(x;\mathrm{\μ},\mathrm{\σ})=\frac{1}{x\sqrt{2\mathrm{\π}{\mathrm{\σ}}^2}}{e}^{-\frac{{(\mathrm{Inx}-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}}---\left(17\right)$

式中μ、σ分别是变量对数的平均值与标准差；

这里采用极大似然估计法计算其未知参数μ和σ：

似然函数为： $L(\mathrm{\μ},\mathrm{\σ})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}{x}_i}{e}^{-\frac{{(\mathrm{In}{x}_i-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}}---\left(18\right)$

两边取对数： $\mathrm{InL}(\mathrm{\μ},{\mathrm{\σ}}^2)=-\frac{n}{2}\mathrm{In}\left(2\mathrm{\π}{\mathrm{\σ}}^2\right)-\mathrm{In}\underset{I=1}{\overset{n}{\mathrm{\Π}}}{x}_i-\frac{1}{2{\mathrm{\σ}}^2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{Inx}}_i-\mathrm{\μ})}^2---\left(19\right)$

似然方程组为 $\left\{\begin{array}{c}\frac{\mathrm{\δInL}(\mathrm{\μ},{\mathrm{\σ}}^2)}{\mathrm{\δ\μ}}=\frac{1}{{\mathrm{\σ}}^2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{Inx}}_i-\mathrm{\μ})}^2=0\\ \frac{\mathrm{\δInL}(\mathrm{\μ},{\mathrm{\σ}}^2)}{\mathrm{\δ\μ}}=-\frac{1}{{2\mathrm{\σ}}^2}+\frac{1}{{2\mathrm{\σ}}^4}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{Inx}}_i-\mathrm{\μ})}^2=0\end{array}\right.---\left(20\right)$

解得： $\mathrm{\μ}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{Inx}}_i,{\mathrm{\σ}}^2=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{Inx}}_i-\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{Inx}}_i)}^2---\left(21\right)$

所述步骤5中，残差平方和指标（Residual Sum of Squares，RSS）的计算公式为：

$\mathrm{RSS}=\mathrm{\Σ}{(f\left(x_i\right)-P_{x_i})}^2---\left(22\right)$

式中，x_i为随机变量历史数据的取值，f(x_i)为x_i对应的拟合函数值，

为变量原始分布中x_i对应的概率值。

通过大量实验数据分析，发现风电功率的不同状态的持续时间特性与逆高斯分布相似，同时，风电场有可能长时间（数小时甚至数日）保持相同的出力状态，且在出力较低或较高的状态保持不变的平均时间要长于保持中等出力状态的平均时间；风电功率持续时间分布特性与逆高斯函数分布曲线大致吻合的这一发现，将有助于对风电功率在某一出力水平下的持续时间做出估计，对于风电功率的短期预测提供重要参考。

本发明的有益效果为：有助于对风电功率在某一出力水平下的持续时间做出估计，对于风电功率的短期预测提供重要参考；对含风电电力系统的规划、运行方式的安排等，都有一定的指导意义。

附图说明

图1为本发明的流程示意图；

图2为Brazos风电场所有状态风电功率持续时间出现次数统计曲线；

图3为Brazos风电场所有状态持续时间分布拟合曲线；

图4为不同风电场风电功率持续时间概率分布特性曲线。

具体实施方式：

下面结合附图与实施例对本发明作进一步说明。

以美国德克萨斯州的Delaware风电场，Brazos风电场和Capridge风电场，澳大利亚的Woolnorth风电场以及爱尔兰、英国与德国TenneT管辖的风电场为例，对风电功率在不同风电场、不同状态下的持续时间特性进行分析、校验；几个风电场的基本信息如表1所示：

表1风电场的基本信息表

实施例步骤如下：

如图1所示：实施步骤（1）：定义风电功率的状态：

如图2所示：以Brazos风电场为例，将风电功率的可能取值范围离散化为11个功率区间，每个功率区间即为风电功率的一个状态；依照定义，实测风电功率序列中的每一个数值均对应一个风电功率状态；该风电场的额定功率为160MW，被等分为10个状态（加0状态共11个状态）。

图中横坐标为状态n的不同持续时间，纵坐标为各持续时间对应的次数，可以看出，每个状态持续时间分布的趋势大致相同。

实施步骤（2）：进行风电功率状态持续时间的概率统计：

风电功率状态持续时间的统计包括两个方面：保持某个状态的时间长度和持续某一个时间长度的次数；例如：当风电功率从任意状态m，进入到状态n后，其中m≠n，开始记录风电功率保持在状态n内的时间；若风电功率经历时间T后跳出状态n，则记录状态n持续时间T一次；按照这种方法统计实测风电功率序列，可以得到风电功率在状态n下不同持续时间各自出现的次数，分别除以对应状态n下持续时间的总次数即得到该状态n下持续时间的概率，统计该状态n下所有持续时间T的概率值，则可以得到风电功率在此状态下持续时间的概率分布。

实施步骤（3）：选择适当的分布函数；

选择Birnbaum-Saunders分布、指数分布、逆高斯分布以及对数正态分布作为参考分布函数；

实施步骤（4）：利用最小二乘估计法和极大似然估计法分别对选择的分布函数进行拟合；

如图3所示，可以看出：指数函数的拟合曲线在持续时间较短的区域与原始分布相差较多；Birnbaum-Saunders分布、逆高斯分布和对数正态分布的拟合效果类似。

如图4所示，7座不同风电场的风电功率状态持续时间的概率分布拟合结果，在此图的状态持续时间的概率分布拟合过程中，每座风电场的出力范围被等分为5个状态，图中只给出其中第1、3、4个状态的分布情况。

对比图3和图4可以看出，当Brazos风电场的有功功率被分为5个状态时（不含0状态），其状态持续时间的概率分布特性与10状态的持续时间概率分布类似。

这说明该风电场的状态持续时间概率分布不受状态定义的影响，均满足概率值随持续时间的增加而逐渐降低的特点，这也表明了逆高斯分布用于拟合状态分布特性具有普适性。

其他风电场状态持续时间的概率分布都具有与Brazos风电场类似的特性，但随着风电场群分布范围越来越广时，由于地域间的相关性影响，其状态持续特性虽然依然大致服从逆高斯分布，但其规律性并没有单一风电场强；例如德国TenneT辖区风电场群的有功功率在0.6p.u.到0.8p.u.之间持续时间的概率分布并不严格满足随持续时间的增加而递减的规律。

因此，上述状态持续时间的概率分布特性，主要适用于单个风电场或小范围内的风电场群。

实施步骤（5）：采用残差平方和指标计算所选择函数的拟合精度；定量比较不同函数拟合效果，选择最优拟合函数。

采用RSS计算上述各函数的拟合精度残差平方和指标（Residual Sum of Squares，RSS）的计算公式为：

$\mathrm{RSS}=\mathrm{\Σ}{(f\left(x_i\right)-P_{x_i})}^2---\left(22\right)$

其中f(x_i)为对应的拟合函数，取为拟合的逆高斯函数；x_i为风电功率序列在某一状态对应持续时间长。

为风电功率序列某一状态持续时间为x_i时实际的概率。

Brazos风电场不同出力状态持续时间的逆高斯分布拟合参数表如表2所示。

其中，

$\mathrm{\μ}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i,\mathrm{\β}={\left(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i\right)}^2\left(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{x_i}\right)-\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i---\left(16\right)$

各函数的拟合精度表如表3所示。

表1Brazos风电场不同出力状态持续时间特性

表3利用四种概率密度函数拟合各状态持续时间的RSS

从表1可以看出，除首末状态外，其余各状态持续时间的λ值均集中在4左右，这进一步验证了上述关于不同状态持续时间分布趋势类似的结论。参数μ随着状态区间功率数值的增加而先减小后增大。这表明风电功率更容易在出力较低或较高的水平保持不变，中等出力水平的平均持续时间较短，可以将其视为风电功率在高、低水平出力之间转换的中间过渡过程。对最长持续时间的统计可以看出，风电功率出力在任意状态下，均有可能持续较长时间。比如在112MW与128MW之间，最长持续时间达到了4259分钟，大约三天的时间。由此可见，虽然风电功率具有波动性和不确定性，但仍然有可能在一个较长的时间段内，保持出力几乎不变。

结合表2、表3可以看出，逆高斯分布拟合的曲线与直方图的RSS最小，因此选用逆高斯分布作为风电功率持续时间分布的拟合函数最为合适。

上述虽然结合附图、具体实施方式对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (3)

1.一种基于极大似然估计的风电功率持续特性描述方法，其特征是：包括以下步骤：

（1）定义风电功率的状态：

将风电功率的可能取值范围离散化为若干个功率区间，每个功率区间即为风电功率的一个状态；依照定义，实测风电功率序列中的每一个数值均对应一个风电功率状态；

（2）进行风电功率状态持续时间的概率统计：

风电功率状态持续时间的统计包括两个方面：保持某个状态的时间长度和持续某一个时间长度的次数；用风电功率在一个状态不同持续时间各自出现的次数，分别除以对应状态持续时间的总次数即得到该状态下持续时间的概率，统计该状态所有持续时间的概率值，则可以得到风电功率在此状态下持续时间的概率分布；

（3）选择适当的分布函数：

选择Birnbaum-Saunders分布、指数分布、逆高斯分布以及对数正态分布作为参考分布函数；

（4）利用最小二乘估计法和极大似然估计分别对选择的分布函数进行拟合；

（5）采用残差平方和指标计算所选择函数的拟合精度。

2.如权利要求1所述的一种基于极大似然估计的风电功率持续特性描述方法，其特征是：所述步骤1中，风电场的额定装机容量记为P_E，拟划分的状态总数为N，则第n个状态代表的功率区间范围设定为：

$\begin{array}{c}(P_{\mathrm{lower}}^n,P_{\mathrm{upper}}^n]\\ P_{\mathrm{lower}}^n=(n-1)\×\frac{P_E}{N},P_{\mathrm{upper}}^n=n\×\frac{P_E}{N},n=\mathrm{1,2},...,N\end{array}---\left(1\right)$

其中，Pⁿ _upper和Pⁿ _lower分别代表第n个状态所代表的风电功率范围的上下限。

3.如权利要求1所述的一种基于极大似然估计的风电功率持续特性描述方法，其特征是：所述步骤5中，残差平方和指标的计算公式为：

$\mathrm{RSS}=\mathrm{\Σ}{(f\left(x_i\right)-P_{x_i})}^2---\left(22\right)$

式中，x_i为随机变量历史数据的取值，f(x_i)为x_i对应的拟合函数值，

为变量原始分布中x_i对应的概率值。

Images