CN103020462A - 计及复杂尾流效应模型的风电场概率输出功率计算方法 - Google Patents

Abstract

本发明公开了风电场输出功率控制技术领域中的一种计及复杂尾流效应模型的风电场概率输出功率计算方法。包括：初始化参数；建立风电场尾流效应模型；初始化仿真年限n，运用序贯蒙特卡罗仿真方法得到风电场各个风电机组和连接电缆的运行状态曲线；在仿真年限n内，利用威布尔分布和风向玫瑰图分别对风速和风向进行模拟，得到规模均为M的风电场风速和风向序列，同时初始化时刻t；利用各个时刻风电场的风速和风向计算相应时刻风电场非故障风电机组的输出功率，进而计算风电场的总输出功率；运用聚类方法处理风电场的总输出功率序列从而得到风电场概率输出功率曲线。本发明有效提高了风电场输出功率模型的精确度，保证了风电场概率输出模型的准确性。

Description

计及复杂尾流效应模型的风电场概率输出功率计算方法

技术领域

本发明属于风电场输出功率控制技术领域，尤其涉及一种计及复杂尾流效应模型的风电场概率输出功率计算方法。

背景技术

由于国际能源危机、环境问题的加剧恶化，风力发电技术的发展和政府部门的激励，风力发电在过去的几年内得到了飞速的增长。因风速具有随机性、波动性和不可控性，风电场的输出功率波动范围通常较大，且速度较快，使得大规模风电并网会带来电网调峰、无功控制及电压控制的困难，严重影响电网的安全稳定及正常调度。因此，如何把握这类电源的变化规律，建立更为符合风电场实际输出状况的模型成为进一步开发风电所迫切要求解决的问题，也是研究风电场大规模并网对电网影响的基础。同时，在研究并网风电场运行和规划方面的相关问题(如：潮流计算、稳定计算、短路计算和随机生产模拟等)时，也需要确定风电场的输出功率。

目前，在建立风电场功率输出模型时，大多数研究人员只考虑简单的尾流效应模型；在建立风电场可靠性模型时，只涉及风电机组的故障率。

发明内容

本发明的目的在于，提供一种计及复杂尾流效应模型的风电场概率输出功率计算方法，用于解决现有技术存在的不足。

为了实现上述目的，本发明提出的技术方案是，一种计及复杂尾流效应模型的风电场概率输出功率计算方法，其特征是所述方法包括：

步骤1：初始化参数，包括风电机组坐标，风电机组的风轮半径r₀，风电机组的额定功率P_r，风电机组的切入风速v_ci，风电机组的切出风速v_co，风电机组的额定风速v_r，风电场风速的威布尔分布的尺度参数c,风电场风速的威布尔分布的形状参数k,风电场的风向玫瑰图,风电机组和连接电缆的故障率λ和修复时间μ；

步骤2：建立风电场尾流效应模型；

所述模型为： $v_{\mathrm{eq}.j}=v_0(1-\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left(\frac{A_{\mathrm{shad}.\mathrm{ij}}}{A_{\mathrm{rotor}}}\right){(1-\frac{v_j(x_{\mathrm{ij}},h_j)}{v_0})}^2});$

其中，v_eq.j为风电场中下游风电机组WT_j在所有上游风电机组共同影响下的等效风速；

v₀为自然风的风速；

N为影响下游风电机组WT_j的上游风电机组的个数；

A_shadij为上游风电机组WT_i在下游风电机组WT_j处的尾流影响区域的面积与下游风电机组WT_j的风轮扫略面积的重叠部分的面积，且

$A_{\mathrm{shad}.\mathrm{ij}}=r_i^2\left(x_{\mathrm{ij}}\right){\cos }^{-1}\left(\frac{r_i^2\left(x_{\mathrm{ij}}\right)+d_{\mathrm{ij}}^2+\mathrm{\Δ}{h}^2-r_j^2}{2r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)+r_j^2{\cos }^{-1}\left(\frac{r_j^2+d_{\mathrm{ij}}^2+{\mathrm{\Δh}}^2-r_i^2\left(x_{\mathrm{ij}}\right)}{2r_j\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)$

$-r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}\·\sin \left[{\cos }^{-1}\left(\frac{r_i^2\left(x_{\mathrm{ij}}\right)+d_{\mathrm{ij}}^2+{\mathrm{\Δh}}^2-r_j^2}{2r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)\right],$

x_ij为上游风电机组WT_i和下游风电机组WT_j之间的距离；

r_i(x_ij)为上游风电机组WT_i在下游风电机组WT_j处的尾流影响半径，由Hensen模型求得，其计算公式为： $r_i\left(x_{\mathrm{ij}}\right)=\left\{\begin{array}{c}r=r_0+k_{w}x\\ k_w=\frac{1}{2\ln (h/h_0)}\end{array}\right.;$

r为距离风电机组x处的尾流影响半径；

r₀为风电机组的风轮半径；

k_w为尾流下降系数；

h为风电机组风轮所处海拔高度；

h₀为地面粗糙度；

h_i为上游风电机组WT_i风轮所处海拔高度；

h_j为下游风电机组WT_j风轮所处海拔高度；

Δh为上游风电机组WT_i风轮所处海拔高度h_i和下游风电机组WT_j风轮所处海拔高度h_j的差值的绝对值；

r_j为下游风电机组WT_j的风轮半径；

d_ij为上游风机WT_i和下游风机WT_j之间的水平间距；

A_rotor为风电机组的风轮扫略面积，且

v_j(x_ij,h_j)为下游风电机组WT_j在上游风电机组WT_i单独影响下的风速，且

$v_j(x_{\mathrm{ij}},h_j)=v_i[1-(1-\sqrt{1-C_T})\·{\left(\frac{h_i}{h_j}\right)}^{2\mathrm{\α}}\·{\left(\frac{r_j}{r_i\left(x_{\mathrm{ij}}\right)}\right)}^2\·\left(\frac{A_{\mathrm{shad}.\mathrm{ij}}}{A_{\mathrm{rotor}}}\right)]{\left(\frac{h_j}{h_i}\right)}^{\mathrm{\α}};$

v_i为上游风电机组WT_i的风速；

C_T为推力系数；

α为风速随高度变化系数，并且

通过Lissaman尾流效应模型计算，其计算公式为

v(h_i)为在高度h_i处测得的风速，v(h_j)为在高度h_j处测得的风速；

步骤3：初始化仿真年限n，运用序贯蒙特卡罗仿真方法得到风电场各个风电机组和连接电缆的运行状态曲线；

步骤4：在仿真年限n内，利用威布尔分布和风向玫瑰图分别对风速和风向进行模拟，得到规模均为M的风电场风速和风向序列，同时初始化时刻t,M＝n×8760；

步骤5：利用t时刻风电场的风速和风向计算t时刻风电场非故障风电机组的输出功率，进而计算风电场的总输出功率，具体包括步骤5.1-步骤5.4；

步骤5.1：根据风电场的风向，利用坐标变换公式 $\left[\begin{array}{c}{x}_l^{\′}\\ y_l^{\′}\\ z_l^{\′}\end{array}\right]=Q\left[\begin{array}{c}{x}_l\\ y_l\\ z_l\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\δ}& \sin \mathrm{\δ}& 0\\ -\sin \mathrm{\δ}& \cos \mathrm{\δ}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_l\\ y_l\\ z_l\end{array}\right],$ 对t时刻风电场各处的非故障风电机组的坐标进行修正；其中，x_l为风电机组l的经度，y_l为风电机组l的纬度，z_l为风电机组l的高度，δ为t时刻与前一时刻的风向的差值，x′_l为变换后的风电机组l的经度，y′_l为变换后的风电机组l的纬度，z′_l为变换后的风电机组l的高度，Q为变换矩阵；

步骤5.2：利用风电场尾流效应模型计算t时刻风电场各处的非故障风电机组处的风速v；

步骤5.3：根据t时刻风电场各处的非故障风电机组的风速，利用风电机组的输出特性曲线计算t时刻风电场各处的非故障风电机组的输出功率；其中，计算t时刻风电场各处的非故障风电机组的输出功率的计算公式为：

其中， $A=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[v_{\mathrm{ci}}(v_{\mathrm{ci}}+v_r)-{4v}_{\mathrm{ci}}{v}_r{\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3],$

$B=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[4(v_{\mathrm{ci}}+v_r){\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3-({3v}_{\mathrm{ci}}+v_r)],$

$C=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[2-4{\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3];$

步骤5.4：根据t时刻风电场各处连接电缆的运行状态，对与风电场故障电缆相关联的风电机组的输出功率予以删除，进而求得t时刻风电场的总输出功率，并令t＝t+1;

步骤6:判断是否有t＞M，如果t＞M，则执行步骤7；否则，返回步骤5；

步骤7：将M个时刻的风电场的总输出功率按时间先后顺序排列，得到风电场的总输出功率序列，运用聚类方法处理风电场的总输出功率序列从而得到风电场概率输出功率曲线。

所述聚类方法为k-means法。

本发明在建立风电场功率输出模型时，采用复杂地形下尾流效应模型，该模型充分考虑风电场实际地形情况，能更好描述风电机组间的尾流效应；在建立风电场可靠性模型时，考虑了风电机组和连接电缆的故障率，更加符合风电场实际的运行状况；最后，采用聚类方法对风电场概率输出功率进行建模。本发明能有效提高风电场输出功率模型的精确度，从而保证风电场概率输出模型的准确性，为风电场的规划、调度提供相关可靠依据。

附图说明

图1是计及复杂尾流效应模型的风电场概率输出功率计算方法流程图；

图2是风电机组的部分遮挡尾流效应模型图；

图3是风电场风机布局图；

图4是风电机组链形连接方式图；

图5是风电场风向玫瑰图；

图6是风电场风速模拟曲线图；

图7是风电机组和连接电缆的运行状态图；

图8是风电场的实时输出功率曲线图；

图9是风电场的概率输出功率曲线图。

具体实施方式

下面结合附图，对优选实施例作详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。

实施例1

实施例1为本发明的原理说明。图1是计及复杂尾流效应模型的风电场概率输出功率计算方法流程图，如图1所示，本发明提供的方法包括：

步骤1：初始化参数，包括风电机组坐标，风电机组的风轮半径r₀，风电机组的额定功率P_r，风电机组的切入风速v_ci，风电机组的切出风速v_co，风电机组的额定风速v_r，风电场风速的威布尔分布的尺度参数c,风电场风速的威布尔分布的形状参数k,风电场的风向玫瑰图,风电机组和连接电缆的故障率λ和修复时间μ。

步骤2：建立风电场尾流效应模型。

当处在风速下游的风电机组受到上游多台风机影响时，需考虑部分遮挡尾流效应模型。下游风电机组WT_j和上游风机WT_i的布局如图2所示，假设二者相距为x_ij上游风电机组WT_i风轮所处海拔高度为h_i，下游风电机组WT_j风轮所处海拔高度为h_j，则下游风电机组WT_j受到上游风电机组WT_i影响后的风速计算公式如下：

$v_j(x_{\mathrm{ij}},h_j)=v_i[1-(1-\sqrt{1-C_T})\·{\left(\frac{h_i}{h_j}\right)}^{2\mathrm{\α}}\·{\left(\frac{r_j}{r_i\left(x_{\mathrm{ij}}\right)}\right)}^2\·\left(\frac{A_{\mathrm{shad}.\mathrm{ij}}}{A_{\mathrm{rotor}}}\right)]{\left(\frac{h_j}{h_i}\right)}^{\mathrm{\α}}---\left(1\right)$

式中：v_i为上游风电机组WT_i处的风速，C_T为推力系数，α为风速随高度变化系数，一般取α＝1/7，r_j为下游风电机组WT_j的风轮半径，A_rotor为风电机组的风轮扫略面积，且

A_shad.ij为上游风电机组WT_i在下游风电机组WT_j处的尾流影响区域的面积与下游风电机组WT_j的风轮扫略面积的重叠部分的面积（见图2阴影部分），计算公式可由下式求得：

$A_{\mathrm{shad}.\mathrm{ij}}=r_i^2\left(x_{\mathrm{ij}}\right){\cos }^{-1}\left(\frac{r_i^2\left(x_{\mathrm{ij}}\right)+d_{\mathrm{ij}}^2+\mathrm{\Δ}{h}^2-r_j^2}{2r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)+r_j^2{\cos }^{-1}\left(\frac{r_j^2+d_{\mathrm{ij}}^2+{\mathrm{\Δh}}^2-r_i^2\left(x_{\mathrm{ij}}\right)}{2r_j\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)$

$-r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}\·\sin \left[{\cos }^{-1}\left(\frac{r_i^2\left(x_{\mathrm{ij}}\right)+d_{\mathrm{ij}}^2+{\mathrm{\Δh}}^2-r_j^2}{2r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)\right]---\left(2\right)$

式中：Δh为上游风电机组WT_j风轮所处海拔高度h_i和下游风电机组WT_j风轮所处海拔高度h_j的差值的绝对值，即Δh＝|h_i-h_j|；r_i(x_ij)为上游风电机组WT_i在下游风电机组WT_j处的尾流影响半径，由Jensen模型求得，其计算公式为：

$r_i\left(x_{\mathrm{ij}}\right)=\left\{\begin{array}{c}r=r_0+k_{w}x\\ k_w=\frac{1}{2\ln (h/h_0)}\end{array}\right.---\left(3\right)$

式中：r为距离风电机组x处的尾流影响半径，单位：米；r₀为风电机组的风轮半径，单位：米；h为风电机组风轮所处海拔高度，h₀为地面粗糙度，h₀为一常数；k_w为尾流下降系数，其取值根据当地的地形和气候条件进行选取。当风力机接收到的是自然风速时，k_w等于0.04，否则k_w等于0.08。

复杂地形下风电场尾流效应模型中还考虑了Lissaman尾流效应模型，即不同风机安装的地点不同，具有不同的海拔高度，而风速随高度变化而变化，导致风电场的风速分布不均。具体公式如下：

$\frac{v\left(h\right)}{v\left(H_0\right)}={\left(\frac{h}{H_0}\right)}^{\mathrm{\α}}---\left(4\right)$

式中：v(H₀)为在高度为H₀测得的风速，单位：米/秒；v(h)为在高度为h测得的风速，单位：米/秒；α为风速随高度变化系数。依据上述公式（4）的指导，在本发明中，上述公式（1）中的

通过

计算。

在多机组的大型风电场中，处在风向下游的风电机组WT_j很有可能受到上游多个风电机组WT_i的尾流影响，故应考虑风电机组的多重尾流效应模型。模型假设多个上游机组对下游机组的共同影响等于各个上游机组对其单独影响之和，故下游机组WT_j处的等效风速v_eq.j满足下式：

$1-\frac{v_{\mathrm{eq}.j}}{v_0}=\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left(\frac{A_{\mathrm{shad}.\mathrm{ij}}}{A_{\mathrm{rotor}}}\right){(1-\frac{v_j(x_{\mathrm{ij}},h_j)}{v_0})}^2}---\left(5\right)$

式中：v₀为自然风的风速，单位：米/秒；N为影响下游风电机组WT_j的上游风电机组的个数；v_eq.j为风电场中下游风电机组WT_j在所有上游风电机组共同影响下的等效风速，单位：米/秒；v_j(x_ij,h_j)为下游风电机组WT_j在上游风电机组WT_i单独影响下的风速，单位：米/秒。

步骤3：初始化仿真年限n，运用序贯蒙特卡罗仿真方法得到风电场各个风电机组和连接电缆的运行状态曲线。

1、风电机组停运的蒙特卡罗仿真

风电机组停运模型采用运行、停运两状态Mar kov模型，用正常状态持续时间和修复时间来描述。由于风电机组结构相对简单，计划检修时间较短，而且计划检修可以安排在低风速时进行，在可靠性分析时可以不予考虑。一般认为正常运行持续时间和修复时间都服从指数分布，故可采用序贯蒙特卡罗法进行仿真，表达式为：

${\mathrm{\τ}}_1=-\frac{1}{\mathrm{\λ}}\ln {\mathrm{\γ}}_1=-T_{\mathrm{MTTF}}\ln {\mathrm{\γ}}_1$

${\mathrm{\τ}}_2=-\frac{1}{\mathrm{\μ}}\ln {\mathrm{\γ}}_2=-T_{\mathrm{MTTR}}\ln {\mathrm{\γ}}_2$

式中：τ₁和τ₂分别为正常状态持续时间和修复时间；λ为风电机组故障率；μ为风电机组修复率；γ₁和γ₂为[0,1]之间的均匀分布随机数；T_MTTF为平均无故障工作时间；T_MTTR为平均修复时间。

当某一风电机组停运时，不仅该机组的输出功率变为零，同时还对整个风电场的尾流效应有影响，即故障机组相当于不存在，对后面的机组无尾流影响。

2、连接电缆停运的蒙特卡罗仿真

连接电缆停运模型同样采用运行、停运两状态模型。在得知连接电缆故障率和修复率时，亦可采用上述风电机组故障蒙特卡罗仿真计算公式进行仿真。

从陆上风电场和海上风电场的设计经验来看，风电场的内部电气接线方式基本上有以下3种形式：链形、环形、星形。在不同的连接方式，连接电缆故障造成风电机组功率损失是不同的。链形情况下，距离回流母线越近的电缆故障造成的功率损失越大；环形情况下，一条连接电缆故障不会造成功率的损失；星形情况下，每台风力发电机及其电缆故障与否都不影响风电场其它部分的正常运行，并且能够实现独立调节。

步骤4：在仿真年限n内，利用威布尔分布和风向玫瑰图分别对风速和风向进行模拟，得到规模均为M的风电场风速和风向序列，同时初始化时刻t,M＝n×8760。

在本发明中，设定每小时为1个时刻，每天为24个时刻，1年365天共有365×24＝8760个时刻。由于设定的年限为n，所以要先得到M＝n×8760个时刻风电场风速和风向序列，每个序列按照时刻先后顺序排序可以得到风电场风速序列和风向序列。

上述步骤3和步骤4已经是本领域技术人员的常用技术，因此在本发明中不再赘述。

步骤5：利用t时刻风电场的风速和风向计算t时刻风电场非故障风电机组的输出功率，进而计算风电场的总输出功率，具体包括步骤5.1-步骤5.4。

步骤5.1：根据风电场的风向，利用坐标变换公式 $\left[\begin{array}{c}{x}_l^{\′}\\ y_l^{\′}\\ z_l^{\′}\end{array}\right]=Q\left[\begin{array}{c}{x}_l\\ y_l\\ z_l\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\δ}& \sin \mathrm{\δ}& 0\\ -\sin \mathrm{\δ}& \cos \mathrm{\δ}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_l\\ y_l\\ z_l\end{array}\right],$ 对t时刻风电场各处的非故障风电机组的坐标进行修正；其中，x_l为风电机组l的经度，y_l为风电机组l的纬度，z_l为风电机组l的高度，δ为t时刻与前一时刻的风向的差值，x′_l为变换后的风电机组l的经度，y′_l为变换后的风电机组l的纬度，z′_l为变换后的风电机组l的高度，Q为变换矩阵。

步骤5.2：利用风电场尾流效应模型计算t时刻风电场各处的非故障风电机组处的风速v。即利用步骤2中得到的公式（5）计算t时刻风电场各处的非故障风电机组处的风速v。

步骤5.3：根据t时刻风电场各处的非故障风电机组的风速，利用风电机组的输出特性曲线计算t时刻风电场各处的非故障风电机组的输出功率；其中，计算t时刻风电场各处的非故障风电机组的输出功率的计算公式为：

其中， $A=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[v_{\mathrm{ci}}(v_{\mathrm{ci}}+v_r)-{4v}_{\mathrm{ci}}{v}_r{\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3],$

$B=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[4(v_{\mathrm{ci}}+v_r){\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3-({3v}_{\mathrm{ci}}+v_r)],$

$C=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[2-4{\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3].$

P_r、v_ci、v_co和v_r分别为风电机组的额定功率、风电机组的切入风速、风电机组的切出风速和风电机组的额定风速。

步骤5.4：根据t时刻风电场各处连接电缆的运行状态，对与风电场故障电缆相关联的风电机组的输出功率予以删除，进而求得t时刻风电场的总输出功率，并令t＝t+1。

根据风电场连接电缆接线方式，分析风电机组是否因为电缆故障而不能将功率传输至汇流母线处，因电缆故障导致不能将功率传输到汇流母线处的风电机组的输出功率将不被计入风电场总的输出功率。

步骤6:判断是否有t＞M，如果t＞M，则执行步骤7；否则，返回步骤5。

步骤7：将M个时刻的风电场的总输出功率按时间先后顺序排列，得到风电场的总输出功率序列，运用聚类方法处理风电场的总输出功率序列从而得到风电场概率输出功率曲线。

聚类是数据挖掘的一个分支，它是指将物理或者抽象的集合分组成为多个类的过程，每一类又称为一个簇。处在同一个簇中的对象之间具有较高的相似性，不同簇中的对象则相异，k-means法是聚类算法中的一种常用方法。

设共有M个样本X₁,X₂，...,X_M，计第i个样本X_i＝[x_i1,x_i2，...,x_ip]，则第i个样本与第j个样本之间的闵科夫斯基距离定义为：

$d_{\mathrm{ij}}\left(q\right)={\left[\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{|x_{\mathrm{ik}}-x_{\mathrm{jk}}|}^q\right]}^{1/q},i=\mathrm{1,2},...,M;j=\mathrm{1,2},...,M---\left(7\right)$

式中q为正整数。当q＝2时，称之为欧氏距离，即

$d_{\mathrm{ij}}={\left[\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{|x_{\mathrm{ik}}-x_{\mathrm{jk}}|}^2\right]}^{1/2}---\left(8\right)$

样本i的轮廓值为

$S\left(i\right)=\frac{\min \left(b\right)-a}{\max [a,\min \left(b\right)]},i=1,2,...,M---\left(9\right)$

式中：a为样本i与同簇的其他样本之间的平均距离；b为一个向量，其元素是样本i与不同簇的簇内样本之间的平均距离。

轮廓值S(i)的取值范围为[-1,1]，S(i)的值越大，说明样本i的分类越合理，当S(i)＜0时，说明样本i的分类不合理，还有比目前分类更合理的方案。

基于上述k-means法的原理，本发明的步骤7的聚类过程是：

步骤101：先将风电场的总输出功率序列分为k个类，k为设定值。也就是将总输出功率序列中的总输出功率分别分配到k个类中，每个类作为一个样本。

步骤102：对每个样本进行标准化处理，即减去均值，除以标准差。

步骤103：选取k个样本作为初始聚类点。

步骤104：计算每个类的重心，作为聚类中心。

步骤105：分别计算选取的k个样本以外的其他样本与初始聚类点的距离，将它们归入与其距离最近的初始聚类点所属的类。

步骤106：重复上述步骤103-105，直至所有样本不能分配为至。

步骤107：根据上述公式（9）计算S(i)，若S(i)≤0，首先重新选取初始聚类点进行聚类，直至S(i)满足S(i)＞0，若所有的初始聚类点均不能满足，则重新输入k值，进行聚类。

同时，应按时段计算样本数据，当本时段样本数据与上一时段不同时，首先使用上一时段的聚类k值进行聚类，聚类结果若与上一时段相同，则继续使用该聚类结果；若与上一时段的结果不同，则重新输入k值聚类，使用本时段样本数据的聚类结果。

实施例2

实施例2结合实际数据说明本发明的实现过程。如图3所示，风电场含有80台风电机组，风电机组旁的数据是以测风塔为基准的相对风机高度。假设该风电场的内部电气接线方式的形式为链形，10kV汇流母线位于风电场左侧(见图4)，图3中第一排风电机组对应图4中的分别为1#~10#风电机组，相应的电缆亦为1#~10#，以此类推；风电机组的具体参数为：风轮半径r₀=45m，切入风速v_ci=3m/s，切出风速v_co=25m/s，额定风速v_r=12m/s，额定功率P_r＝2MW，轮毂高度h=60m；风电场每月风速的威布尔分布参数c，k如表1所示，风电场的风电玫瑰图如图5所示。根据上述参数数据利用威布尔分布和风向玫瑰图进行模拟得到风电场的风速和风向，仿真时间为3年，仿真结果如图6所示。

表1风电场月风速的威布尔分布表

2、运用序贯蒙特卡罗仿真得到风电场各个风电机组和连接电缆的运行状态曲线，取风电机组故障率λ_WT=0.012次/年，修复率μ_WT=30天，取连接电缆故障率λ_CA=0.008次/年，μ_CA=12天，仿真时间为3年，仿真结果如图7所示，图中只画出了在模拟时间内发生故障的风电机组和连接电缆，未发生故障的没在图中表示。

3、运用复杂地形尾流效应模型和计及电缆故障的风电场可靠性模型，将风电场实时风速和风向代入到上述模型中进行计算，得到风电场的实时输出功率，计算结果如图8所示。

4、运用聚类思想对风电场的输出功率曲线进行聚类，得到风电场的概率输出功率曲线见图9。

本发明是计及复杂尾流效应模型的风电场概率输出功率建模方法，具有以下特点：

（1）采用复杂尾流效应模型模拟风电机组间的尾流影响，更加符合实际情况，有助于准确计算出各个风电机组的输入风速；

（2）考虑风电机组和连接电缆的故障率，贴近风电场的实际运行情况，更能准确描述风电场的整体功率输出；

（3）利用聚类方法建立风电场的概率输出功率模型，能更好地描述风电场的概率输出。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (2)

1.一种计及复杂尾流效应模型的风电场概率输出功率计算方法，其特征是所述方法包括：

步骤1：初始化参数，包括风电机组坐标，风电机组的风轮半径R₀，风电机组的额定功率P_r，风电机组的切入风速v_ci，风电机组的切出风速v_co，风电机组的额定风速v_r，风电场风速的威布尔分布的尺度参数c，风电场风速的威布尔分布的形状参数k，风电场的风向玫瑰图，风电机组和连接电缆的故障率λ和修复时间μ；

步骤2：建立风电场尾流效应模型；

所述模型为： $v_{\mathrm{eq}.j}=v_0(1-\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left(\frac{A_{\mathrm{shad}.\mathrm{ij}}}{A_{\mathrm{rotor}}}\right){(1-\frac{v_j(x_{\mathrm{ij}},h_j)}{v_0})}^2});$

其中，v_eq.j为风电场中下游风电机组WT_j在所有上游风电机组共同影响下的等效风速；

v₀为自然风的风速；

N为影响下游风电机组WT_j的上游风电机组的个数；

A_shad.ij为上游风电机组WT_i在下游风电机组WT_j处的尾流影响区域的面积与下游风电机组WT_j的风轮扫略面积的重叠部分的面积，且

$A_{\mathrm{shad}.\mathrm{ij}}=r_i^2\left(x_{\mathrm{ij}}\right){\cos }^{-1}\left(\frac{r_i^2\left(x_{\mathrm{ij}}\right)+d_{\mathrm{ij}}^2+\mathrm{\Δ}{h}^2-r_j^2}{2r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)+r_j^2{\cos }^{-1}\left(\frac{r_j^2+d_{\mathrm{ij}}^2+{\mathrm{\Δh}}^2-r_i^2\left(x_{\mathrm{ij}}\right)}{2r_j\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)$

$-r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}\·\sin \left[{\cos }^{-1}\left(\frac{r_i^2\left(x_{\mathrm{ij}}\right)+d_{\mathrm{ij}}^2+{\mathrm{\Δh}}^2-r_j^2}{2r_i\left(x_{\mathrm{ij}}\right)\·\sqrt{{\mathrm{\Δh}}^2+d_{\mathrm{ij}}^2}}\right)\right],$

x_ij为上游风电机组WT_i和下游风电机组WT_j之间的距离；

r_i(x_ij)为上游风电机组WT_i在下游风电机组WT_j处的尾流影响半径，由Jensen模型求得，其计算公式为： $r_i\left(x_{\mathrm{ij}}\right)=\left\{\begin{array}{c}r=r_0+k_{w}x\\ k_w=\frac{1}{2\ln (h/h_0)}\end{array}\right.;$

r为距离风电机组x处的尾流影响半径；

r₀为风电机组的风轮半径；

k_w为尾流下降系数；

h为风电机组风轮所处海拔高度；

h₀为地面粗糙度；

h_i为上游风电机组WT_i风轮所处海拔高度；

h_j为下游风电机组WT_j风轮所处海拔高度；

Δh为上游风电机组WT_i风轮所处海拔高度h_i和下游风电机组WT_j风轮所处海拔高度h_j的差值的绝对值；

r_j为下游风电机组WT_j的风轮半径；

d_ij为上游风机WT_i和下游风机WT_j之间的水平间距；

A_rotor为风电机组的风轮扫略面积，且

v_j(x_ij,h_j)为下游风电机组WT_j在上游风电机组WT_i单独影响下的风速，且 $v_j(x_{\mathrm{ij}},h_j)=v_i[1-(1-\sqrt{1-C_T})\·{\left(\frac{h_i}{h_j}\right)}^{2\mathrm{\α}}\·{\left(\frac{r_j}{r_i\left(x_{\mathrm{ij}}\right)}\right)}^2\·\left(\frac{A_{\mathrm{shad}.\mathrm{ij}}}{A_{\mathrm{rotor}}}\right)]{\left(\frac{h_j}{h_i}\right)}^{\mathrm{\α}};$

v_i为上游风电机组WT_i的风速；

C_T为推力系数；

α为风速随高度变化系数，并且

通过Lissaman尾流效应模型计算，其计算公式为

v(h_i)为在高度h_i处测得的风速，v(h_j)为在高度h_j处测得的风速；

步骤3：初始化仿真年限n，运用序贯蒙特卡罗仿真方法得到风电场各个风电机组和连接电缆的运行状态曲线；

步骤4：在仿真年限n内，利用威布尔分布和风向玫瑰图分别对风速和风向进行模拟，得到规模均为M的风电场风速和风向序列，同时初始化时刻t,M＝n×8760；

步骤5：利用t时刻风电场的风速和风向计算t时刻风电场非故障风电机组的输出功率，进而计算风电场的总输出功率，具体包括步骤5.1-步骤5.4；

步骤5.1：根据风电场的风向，利用坐标变换公式 $\left[\begin{array}{c}{x}_l^{\′}\\ y_l^{\′}\\ z_l^{\′}\end{array}\right]=Q\left[\begin{array}{c}{x}_l\\ y_l\\ z_l\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\δ}& \sin \mathrm{\δ}& 0\\ -\sin \mathrm{\δ}& \cos \mathrm{\δ}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_l\\ y_l\\ z_l\end{array}\right],$ 对t时刻风电场各处的非故障风电机组的坐标进行修正；其中，x_l为风电机组l的经度，y_l为风电机组l的纬度，z_l为风电机组l的高度，δ为t时刻与前一时刻的风向的差值，x′_l为变换后的风电机组l的经度，y′_l为变换后的风电机组l的纬度，z′_l为变换后的风电机组l的高度，Q为变换矩阵；

步骤5.2：利用风电场尾流效应模型计算t时刻风电场各处的非故障风电机组处的风速v；

步骤5.3：根据t时刻风电场各处的非故障风电机组的风速，利用风电机组的输出特性曲线计算t时刻风电场各处的非故障风电机组的输出功率；其中，计算t时刻风电场各处的非故障风电机组的输出功率的计算公式为：

其中， $A=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[v_{\mathrm{ci}}(v_{\mathrm{ci}}+v_r)-{4v}_{\mathrm{ci}}{v}_r{\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3],$

$B=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[4(v_{\mathrm{ci}}+v_r){\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3-({3v}_{\mathrm{ci}}+v_r)],$

$C=\frac{1}{{(v_{\mathrm{ci}}-v_r)}^2}[2-4{\left(\frac{v_{\mathrm{ci}}+v_r}{{2v}_r}\right)}^3];$

步骤5.4：根据t时刻风电场各处连接电缆的运行状态，对与风电场故障电缆相关联的风电机组的输出功率予以删除，进而求得t时刻风电场的总输出功率，并令t＝t+1;

步骤6:判断是否有t＞M，如果t＞M，则执行步骤7；否则，返回步骤5；

步骤7：将M个时刻的风电场的总输出功率按时间先后顺序排列，得到风电场的总输出功率序列，运用聚类方法处理风电场的总输出功率序列从而得到风电场概率输出功率曲线。

2.根据权利要求1所述的方法，其特征是所述聚类方法为k-means法。

Images