CN103208798A - 一种含风电场电力系统概率潮流的计算方法 - Google Patents

Abstract

本发明公开了一种含风电场电力系统概率潮流的计算方法，包括如下步骤：步骤1：计算风电场功率以及负荷功率的概率分布；步骤2：用牛顿法进行含风电场电力系统确定性潮流的计算，求出灵敏度矩阵S₀；步骤3：计算各节点注入向量的各阶半不变量；步骤4：根据注入量的各阶半不变量分别求解状态变量的各阶半不变量；步骤5：根据Gram-Charlier级数展开求解概率密度函数和累积分布函数。本发明可以有效减少常规的将风电场处理为功率因数恒定的简化PQ模型时半不变量法求解带来的误差。

Description

一种含风电场电力系统概率潮流的计算方法

技术领域

本发明涉及一种含风电场电力系统概率潮流的计算方法，具体涉及一种基于异步风力发电机（简称“异步风机”或“异步机”）RX模型的含风电场电力系统概率潮流的计算方法，适用于含异步风力发电机的含风电场电力系统。本发明可以减小常规的将风电场处理为功率因数恒定的简化PQ模型时半不变量法求解的误差。

背景技术

随着世界范围内大规模风电并网运行，风电并网对电网的影响研究成为当前电力系统的研究热点。风能的随机性使得确定性分析方法提供的结果过于保守，计及随机性进行风电并网研究可以为调度提供更可靠的参考，因此运用概率技术分析风电并网问题具有十分重要的意义。

含风电场电力系统概率潮流计算能够计及风能的随机性对并网风电系统进行评估。文献一《考虑尾流效应的风电场建模拟及随机潮流计算》（西安交通大学学报，2008年第42卷第12期第1515页）应用蒙特卡罗法对含风电场系统进行了概率潮流计算，分析了风电场并网运行对系统各节点电压的影响，但该方法需要多次模拟，计算时间较长。文献二《计及分布式发电的配电系统随机潮流计算》（电力系统自动化2005年第29卷第24期第15页）重点研究了分布式发电中的风力发电和太阳能发电的随机出力对配电系统电压的影响，文献三《含风电场的电力系统概率潮流计算》（电网技术2009年第33卷第16期第87页）通过建立风电机组概率模型，运用概率潮流计算分析了风电场加入前后系统电压与潮流的变化情况。二者都采用基于半不变量法和Gram-Charlier级数展开法获得状态量的概率分布，计算效率比蒙特卡罗法大大提高。文献四《基于半不变量法的随机潮流误差分析》（电网技术2009年第33卷第18期第32页）指出运用半不变量法进行随机潮流计算时，通常需要做节点注入功率之间相互独立的假设，对于节点注入功率线性相关的情况，需要做特殊处理，否则会引起较大误差。而文献二与文献三在处理风电场节点的有功和无功时，假设功率因数恒定，即Q＝Ptanα，也即风电场节点注入的有功和无功之间线性相关，可能会带来较大误差。

因此，传统的将风电场处理为功率因数恒定的简化PQ模型可能带来较大误差。

发明内容

发明目的：针对上述现有技术存在的将风电场处理为功率因数恒定的简化PQ模型时不满足半不变量法要求的节点注入量之间相互独立的条件从而带来的误差较大的问题和不足，本发明的目的是提供一种含风电场电力系统概率潮流的计算方法，可以有效减少常规的将风电场处理为功率因数恒定的简化PQ模型时半不变量法求解带来的误差。

技术方案：为实现上述发明目的，本发明采用的技术方案为一种含风电场电力系统概率潮流的计算方法，包括如下步骤：

步骤1：计算风电场功率以及负荷功率的概率分布；

步骤2：用牛顿法进行含风电场电力系统确定性潮流的计算，求出灵敏度矩阵S₀；

步骤3：计算各节点注入向量的各阶半不变量；

步骤4：根据注入量的各阶半不变量分别求解状态变量的各阶半不变量；

步骤5：根据Gram-Charlier级数展开求解概率密度函数和累积分布函数。

进一步的，所述步骤1包括：

假设风速服从三参数的威布尔分布，则风速的分布函数F_Weibull(v)为：

$F_{\mathrm{Weibull}}\left(v\right)=1-\exp [-{\left(\frac{v-v_0}{c}\right)}^k]$

式中，v为风速，k、c、v₀分别代表威布尔分布的三个参数，其中k代表形状参数，c代表尺度参数，v₀代表位置参数；

根据给定风机的切入、切出及额定风速，计算风机的输出功率P_W如下式：

$P_W=\left\{\begin{array}{cc}0,& v\&le;v_{\mathrm{ci}}\\ k_{1}v+k_2,& v_{\mathrm{ci}}<v{\&le;v}_r\\ P_r,& v_r<v\&le;v_{\mathrm{co}}\\ 0,& v>v_{\mathrm{co}}\end{array}\right.$

式中：k₁＝P_r/(v_r-v_ci)；k₂＝-k₁v_ci，P_r为风力发电机的额定功率，v_ci为切入风速，v_r为额定风速，v_co为切出风速；

分别根据负荷功率数据x₁、x₂、…、x_n和风机的输出功率P_W的历史数据，计算期望和方差。

进一步的，所述步骤2包括：

假设节点i连接有风电场，则与节点i对应的潮流方程为：

$\left\{\begin{array}{c}{P}_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)-P_{\mathrm{Li}}-V_i\underset{j\&Element;i}{\mathrm{\&Sigma;}}{V}_j(G_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})=0\\ Q_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)-Q_{\mathrm{Li}}-V_i\underset{j\&Element;i}{\mathrm{\&Sigma;}}{V}_j(G_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}})=0\\ P_{\mathrm{mi}}-P_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)=0\end{array}\right.$

式中，P_ei(V_i,θ_i,s_i)、Q_ei(V_i,θ_i,s_i)分别表示与变量V_i、θ_i和s_i有关的风机有功功率和无功功率；P_Li、Q_Li分别表示负荷有功功率和无功功率；V_i、θ_i分别为节点i的电压幅值和相角；G_ij、B_ij分别为节点i、j间的电导与电纳值；θ_ij为节点i、j间的相角差；s_i为异步机滑差；P_mi为风力发电机的机械功率。

求解上述潮流方程：在已知风机有功功率和无功功率P_ei(V_i,θ_i,s_i)、Q_ei(V_i,θ_i,s_i)，负荷有功功率和无功功率P_Li、Q_Li，异步机滑差s_i，G_ij、B_ij的情况下，根据上述潮流方程，计算风电场节点电压幅值V_i和相角θ_i，通用表达式如下：

W＝f(X)

式中，W为节点注入向量，包括P_ei(V_i,θ_i,s_i)和Q_ei(V_i,θ_i,s_i)，为已知；X为节点状态变量，包括节点的电压幅值V_i和相角θ_i，为未知；

在概率潮流计算中，节点注入向量是随机变量，因此将其表示为：

W＝W₀+ΔW

式中，W₀为节点注入向量W的期望值，ΔW为节点注入向量W的随机扰动；

同理将状态变量写成

X＝X₀+ΔX

式中，X₀是状态变量X的期望值，ΔX为状态变量X的随机变量；

进行泰勒级数展开，忽略高次项，得到：

W＝W₀+ΔW＝f(X₀+ΔX)＝f(X₀)+J₀ΔX

其中：

W₀＝f(X₀)

通过牛顿法潮流计算解出系统的状态变量X₀，J₀为牛顿法潮流计算最后一次迭代使用的雅克比矩阵；

ΔW＝J₀ΔX

因此有：

$\mathrm{\&Delta;X}={J_0}^{-1}\mathrm{\&Delta;W}=S_0\mathrm{\&Delta;W}$

式中，S₀为J₀的逆矩阵，称为灵敏度矩阵。

进一步的，所述步骤3包括：

分别计算负荷功率数据x₁、x₂、…、x_n的k阶中心矩，并记作β_k，计算方法如下：

${\mathrm{\&beta;}}_k=\underset{i=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(x_i-m_k)}^k{p}_i$

式中，m_x为负荷功率的期望，p_i为x₁、x₂、…、x_n的分布概率，其中i=1、2、……、n，各阶半不变量由各自的各阶中心矩算得，通用计算方式如下：

γ₁＝m_x

${\mathrm{\&gamma;}}_2={{\mathrm{\&delta;}}_x}^2$

γ₃＝β₃

${\mathrm{\&gamma;}}_4={\mathrm{\&beta;}}_4-3{\mathrm{\&beta;}}_2^2$

γ₅＝β₅-10β₂β₃

${\mathrm{\&gamma;}}_6={\mathrm{\&beta;}}_6-15{\mathrm{\&beta;}}_2{\mathrm{\&beta;}}_4-10{\mathrm{\&beta;}}_3^2+30{\mathrm{\&beta;}}_2^3$

式中，γ_i代表i阶半不变量，δ_x为负荷功率的方差；

根据上述公式，由各变量的已知数据，分别计算节点i注入有功功率和无功功率的n阶半不变量

和

风电场节点的注入有功功率和无功功率的n阶半不变量

和

分别计算滑差与电磁功率不平衡量的n阶半不变量和机械功率与电磁功率不平衡量的n阶半不变量。

进一步的，所述步骤4包括：

根据下式，由注入量的各阶半不变量求得状态量的半不变量：

$\left[\begin{array}{c}{\mathrm{\&Delta;\&theta;}}_1^{\left(n\right)}\\ {\mathrm{\&Delta;v}}_1^{\left(n\right)}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;\&theta;}}_n^{\left(n\right)}\\ {\mathrm{\&Delta;v}}_n^{\left(n\right)}\\ {\mathrm{\&Delta;s}}_n^{\left(n\right)}\end{array}\right]=\left[\begin{array}{ccccc}{S}_{\mathrm{1,1}}^n& S_{\mathrm{1,2}}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{\mathrm{1,2}n}^n& S_{\mathrm{1,2}n+1}^n\\ S_{\mathrm{2,1}}^n& S_{\mathrm{2,2}}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{\mathrm{2,2}n}^n& S_{\mathrm{2,2}n+1}^n\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;\\ S_{2n-\mathrm{1,1}}^n& S_{2n-\mathrm{1,2}}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n-\mathrm{1,2}n}^n& S_{2n-\mathrm{1,2}n+1}^n\\ S_{2n,1}^n& S_{2n,2}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n,2n}^n& S_{2n,2n+1}^n\\ S_{2n+\mathrm{1,1}}^n& S_{2n+1,2}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n+\mathrm{1,2}n}^n& S_{2n+\mathrm{1,2}n+1}^n\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Delta;P}}_1^{\left(n\right)}\\ {\mathrm{\&Delta;Q}}_1^{\left(n\right)}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}\\ {\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}\\ {\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}\end{array}\right]$

式中，

和

表示节点i的电压幅值和相角的n阶半不变量，

和

分别表示滑差与电磁功率不平衡量的n阶半不变量和机械功率与电磁功率不平衡量的n阶半不变量；

表示灵敏度矩阵第i行第j列元素的n次方；

通过上式推导，状态变量的半不变量通过下式计算得到：

${\mathrm{\&Delta;\&theta;}}_i^{\left(n\right)}=S_{2i-\mathrm{1,1}}^n{\mathrm{\&Delta;P}}_1^{\left(n\right)}+S_{2i-\mathrm{1,2}}^n{\mathrm{\&Delta;Q}}_1^{\left(n\right)}+\&CenterDot;\&CenterDot;\&CenterDot;+\&CenterDot;\&CenterDot;\&CenterDot;S_{2i-\mathrm{1,2}n-1}^n{\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}+S_{2i-\mathrm{1,2}n}^n{\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}+S_{2i-\mathrm{1,2}n+1}^n{\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}$

${\mathrm{\&Delta;v}}_i^{\left(n\right)}=S_{2i,1}^n{\mathrm{\&Delta;P}}_1^{\left(n\right)}+S_{2i,2}^n{\mathrm{\&Delta;Q}}_1^{\left(n\right)}+\&CenterDot;\&CenterDot;\&CenterDot;+\&CenterDot;\&CenterDot;\&CenterDot;S_{2i,2n-1}^n{\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}+S_{2i,2n}^n{\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}+S_{2i,2n+1}^n{\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}$

进一步的，所述步骤5包括：

通过节点i的电压幅值和相角的n阶半不变量和反推其各阶中心距，其统一表达式为：

γ₁＝m_x

${\mathrm{\&gamma;}}_2={{\mathrm{\&delta;}}_x}^2$

γ₃＝β₃

${\mathrm{\&gamma;}}_4={\mathrm{\&beta;}}_4-3{\mathrm{\&beta;}}_2^2$

γ₅＝β₅-10β₂β₃

${\mathrm{\&gamma;}}_6={\mathrm{\&beta;}}_6-15{\mathrm{\&beta;}}_2{\mathrm{\&beta;}}_4-10{\mathrm{\&beta;}}_3^2+30{\mathrm{\&beta;}}_2^3$

分布函数F(ξ)的统一表达式为:

F(ξ)＝c₀Φ(ξ)+c₁Φ′(ξ)+c₂Φ′′(ξ)+c₃Φ⁽³⁾(ξ)+c₄Φ⁽⁴⁾(ξ)+c₅Φ⁽⁵⁾(ξ)+c₆Φ⁽⁶⁾(ξ)式中，Φ(ξ)为标准正态分布函数表达式，各系数由下式求得

c₀＝1

c₁＝c₂＝0

$c_3=\frac{1}{3!}(-\frac{{\mathrm{\&beta;}}_3}{{\mathrm{\&delta;}}_3})$

$c_3=\frac{1}{3!}(-\frac{{\mathrm{\&beta;}}_3}{{\mathrm{\&delta;}}_3})$

$c_4=\frac{1}{4!}(-\frac{{\mathrm{\&beta;}}_4}{{\mathrm{\&delta;}}_4})$

$c_5=\frac{1}{5!}(-\frac{{\mathrm{\&beta;}}_5}{{\mathrm{\&delta;}}_5})$

$c_6=\frac{1}{6!}(-\frac{{\mathrm{\&beta;}}_6}{{\mathrm{\&delta;}}_6})$

将电压幅值和相角的各阶中心距代入，即求解出电压幅值和相角的分布函数。

有益效果：含风电场电力系统的概率潮流计算是当前的热点之一。一般假设风电场功率因数恒定，即无功与有功之间具有线性关系，这不满足半不变量法要求的节点注入随机变量之间相互独立的假设，会带来较大误差。本发明采用考虑滑差的异步机RX模型，将风电场有功和无功描述成电压和滑差的函数，采用基于半不变量和Gram-Charlier级数展开法求解状态变量的概率分布，能够有效消除上述误差，进一步提高计算精度。

附图说明

图1是风速的威布尔分布图；

图2是风功率曲线图；

图3是计及滑差的异步发电机等效电路与功率传递关系图；

图4是计及滑差的异步发电机简化等值电路图；

图5是IEEE-14节点系统加入12MW风电场后14号节点的电压分布函数图；

图6是IEEE-14节点系统加入18MW风电场后14号节点的电压分布函数图；

图7是IEEE-14节点系统加入24MW风电场后14号节点的电压分布函数图；

图8是IEEE-30节点系统加入12MW风电场后29号节点的电压分布函数图；

图9是IEEE-30节点系统加入18MW风电场后29号节点的电压分布函数图；

图10是IEEE-30节点系统加入24MW风电场后29号节点的电压分布函数图。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

本发明提出的基于异步机RX模型的含风电场电力系统概率潮流计算，包括以下步骤：

步骤1：计算风电场功率以及负荷功率的概率分布。

步骤2：用牛顿法进行含风电场电力系统确定性潮流计算，求出灵敏度矩阵S₀。

步骤3：计算各节点注入量的各阶半不变量。

步骤4：根据注入量的各阶半不变量分别求解状态变量（简称“状态量”）的各阶半不变量。

步骤5：根据Gram-Charlier级数展开求解概率密度函数和累积分布函数。

所述步骤1提供系统的各项分布数据，具体方法为：

假设风速服从三参数的威布尔分布，则风速的分布函数为：

$F_{\mathrm{Weibull}}\left(v\right)=1-\exp [-{\left(\frac{v-v_0}{c}\right)}^k]$

式中，v为风速，k、c、v₀分别代表威布尔分布的三个参数。其中k代表形状参数，它反映了风速的分布特点，c代表尺度参数，反映的是该地区的平均风速大小，v₀代表位置参数。

根据给定风机的切入、切出及额定风速，计算风机的输出功率如下式：

$P_W=\left\{\begin{array}{cc}0,& v\&le;v_{\mathrm{ci}}\\ k_{1}v+k_2,& v_{\mathrm{ci}}<v{\&le;v}_r\\ P_r,& v_r<v\&le;v_{\mathrm{co}}\\ 0,& v>v_{\mathrm{co}}\end{array}\right.$

式中：k₁＝P_r/(v_r-v_ci)；k₂＝-k₁v_ci；P_r为风力发电机的额定功率；v_ci为切入风速；v_r为额定风速；v_co为切出风速。

上述数据一般由风电场直接给出。

根据负荷功率的历史数据，计算期望和方差，通用公式如下：

如负荷功率数据x₁…x_n，对于服从正态分布的负荷，各数值的分布概率为p₁…p_n，计算负荷的期望m_x和方差δ_x：

$m_x=\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i\right)/n$

${\mathrm{\&delta;}}_x=\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_i-m_x)}^2\right)/n$

同理，将负荷风电场输出功率数据P_W代入，求解相应的期望与方差。

所述步骤2为基于RX模型的含风电场系统潮流计算，具体方法为：

假设节点i连接有风电场，则与节点i对应的潮流方程为：

$\left\{\begin{array}{c}{P}_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)-P_{\mathrm{Li}}-V_i\underset{j\&Element;i}{\mathrm{\&Sigma;}}{V}_j(G_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})=0\\ Q_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)-Q_{\mathrm{Li}}-V_i\underset{j\&Element;i}{\mathrm{\&Sigma;}}{V}_j(G_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}})=0\\ P_{\mathrm{mi}}-P_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)=0\end{array}\right.$

式中，P_ei、Q_ei分别表示风机有功功率和无功功率；P_Li、Q_Li分别表示负荷有功功率和无功功率；V_i、θ_i为节点i的电压幅值和相角；G_ij、B_ij分别为节点i、j间的电导与电纳值；θ_ij为节点i、j间的相角差；s_i为异步机滑差；P_mi为风力发电机的机械功率。

在上述潮流方程中，将未知量V_i、θ_i移至等式右侧，进行方程求解。

即，在已知风机有功功率和无功功率P_ei、Q_ei，负荷有功功率和无功功率P_Li、Q_Li，异步机滑差s_i，节点导纳G_ij、B_ij的情况下，根据上述潮流方程，计算风电场节点电压幅值和相角V_i、θ_i，通用表达式如下：

W＝f(X)

式中，W为节点注入向量，包括节点注入有功功率及无功功率P_ei、Q_ei，为已知。X为节点状态变量，包括节点的电压幅值和相角V_i、θ_i，为未知。

在概率潮流计算中，节点注入量是随机变量，因此可以将其表示为：

W＝W₀+ΔW

式中，W₀为节点注入向量W的期望值，ΔW为节点注入向量W的随机扰动。

同理可将状态变量写成

X＝X₀+ΔX

式中，X₀是状态变量X的期望值，ΔX为状态变量X的的随机变量。

进行泰勒级数展开，忽略高次项，得到：

W＝W₀+ΔW＝f(X₀+ΔX)＝f(X₀)+J₀ΔX

其中：

W₀＝f(X₀)

可以通过牛顿法潮流计算解出系统的状态变量X₀，J₀为牛顿法潮流计算最后一次迭代使用的雅克比矩阵。

ΔW＝J₀ΔX

因此有：

$\mathrm{\&Delta;X}={J_0}^{-1}\mathrm{\&Delta;W}=S_0\mathrm{\&Delta;W}$

式中，S₀为J₀的逆矩阵，称为灵敏度矩阵。

通过上述变换，计算得到灵敏度矩阵S₀各元素。

潮流计算及灵敏度矩阵的详细表述方法见《现代电力系统分析》（2003年，科学出版社）。

所述步骤3为各计算节点注入量的各阶半不变量，统一计算方法具体为：

如负荷功率数据x₁…x_n，计算其k阶中心矩，并记作β_k，计算方法如下：

${\mathrm{\&beta;}}_k=\underset{i=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(x_i-m_k)}^k{p}_i$

各阶半不变量可以由各自的各阶中心矩算得，通用计算方式如下：

γ1_＝m_x

${\mathrm{\&gamma;}}_2={{\mathrm{\&delta;}}_x}^2$

γ₃＝β₃

${\mathrm{\&gamma;}}_4={\mathrm{\&beta;}}_4-3{\mathrm{\&beta;}}_2^2$

γ₅＝β₅-10β₂β₃

${\mathrm{\&gamma;}}_6={\mathrm{\&beta;}}_6-15{\mathrm{\&beta;}}_2{\mathrm{\&beta;}}_4-10{\mathrm{\&beta;}}_3^2+30{\mathrm{\&beta;}}_2^3$

式中，γ_i代表i阶半不变量。

根据上述统一公式，由各变量的已知数据，计算节点i注入有功功率和无功功率的n阶半不变量

风电场节点的注入有功功率和无功功率的n阶半不变量

分别计算滑差与电磁功率不平衡量的n阶半不变量和机械功率与电磁功率不平衡量的n阶半不变量。

所述步骤4为由节点注入功率的各阶半不变量分别求取状态变量的各阶半不变量，具体步骤为：

根据下式，由注入量的各阶半不变量求得状态量

的半不变量：

$\left[\begin{array}{c}{\mathrm{\&Delta;\&theta;}}_1^{\left(n\right)}\\ {\mathrm{\&Delta;v}}_1^{\left(n\right)}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;\&theta;}}_n^{\left(n\right)}\\ {\mathrm{\&Delta;v}}_n^{\left(n\right)}\\ {\mathrm{\&Delta;s}}_n^{\left(n\right)}\end{array}\right]=\left[\begin{array}{ccccc}{S}_{11}^n& S_{12}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{\mathrm{1,2}n}^n& S_{\mathrm{1,2}n+1}^n\\ S_{21}^n& S_{22}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{\mathrm{2,2}n}^n& S_{\mathrm{2,2}n+1}^n\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;\\ S_{2n-\mathrm{1,1}}^n& S_{2n-\mathrm{1,2}}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n-\mathrm{1,2}n}^n& S_{2n-\mathrm{1,2}n+1}^n\\ S_{2n,1}^n& S_{2n,2}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n,2n}^n& S_{2n,2n+1}^n\\ S_{2n+\mathrm{1,1}}^n& S_{2n+1,2}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n+\mathrm{1,2}n}^n& S_{2n+\mathrm{1,2}n+1}^n\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Delta;P}}_1^{\left(n\right)}\\ {\mathrm{\&Delta;Q}}_1^{\left(n\right)}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}\\ {\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}\\ {\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}\end{array}\right]$

式中，

表示节点i的电压幅值和相角的n阶半不变量，

表示滑差及机械功率与电磁功率不平衡量的n阶半不变量；

表示灵敏度矩阵第i行第j列元素的n次方。

通过上式推导可得，状态量的半不变量可通过下式代数表达式计算得到：

${\mathrm{\&Delta;\&theta;}}_k^{\left(n\right)}=S_{2k-\mathrm{1,1}}^n{\mathrm{\&Delta;P}}_1^{\left(n\right)}+S_{2k-\mathrm{1,2}}^n{\mathrm{\&Delta;Q}}_1^{\left(n\right)}+\&CenterDot;\&CenterDot;\&CenterDot;+\&CenterDot;\&CenterDot;\&CenterDot;S_{2k-\mathrm{1,2}n-1}^n{\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}+S_{2k-\mathrm{1,2}n}^n{\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}+S_{2k-\mathrm{1,2}n+1}^n{\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}$

${\mathrm{\&Delta;v}}_k^{\left(n\right)}=S_{2k,1}^n{\mathrm{\&Delta;P}}_1^{\left(n\right)}+S_{2k,2}^n{\mathrm{\&Delta;Q}}_1^{\left(n\right)}+\&CenterDot;\&CenterDot;\&CenterDot;+\&CenterDot;\&CenterDot;\&CenterDot;S_{2k,2n-1}^n{\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}+S_{2k,2n}^n{\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}+S_{2k,2n+1}^n{\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}$

所述步骤5根据Gram-Charlier级数展开求得状态量的分布函数。具体步骤为：

通过节点i的电压幅值和相角的n阶半不变量

反推其各阶中心距，其统一表达式为：

γ₁＝m_x

${\mathrm{\&gamma;}}_2={{\mathrm{\&delta;}}_x}^2$

γ₃＝β₃

${\mathrm{\&gamma;}}_4={\mathrm{\&beta;}}_4-3{\mathrm{\&beta;}}_2^2$

γ₅＝β₅-10β₂β₃

${\mathrm{\&gamma;}}_6={\mathrm{\&beta;}}_6-15{\mathrm{\&beta;}}_2{\mathrm{\&beta;}}_4-10{\mathrm{\&beta;}}_3^2+30{\mathrm{\&beta;}}_2^3$

分布函数的统一表达式为:

F(ξ)＝c₀Φ(ξ)+c₁Φ′(ξ)+c₂Φ′′(ξ)+c₃Φ⁽³⁾(ξ)+c₄Φ⁽⁴⁾(ξ)+c₅Φ⁽⁵⁾(ξ)+c₆Φ⁽⁶⁾(ξ)式中，Φ(ξ)为标准正态分布函数表达式，各系数由下式求得

c₀＝1

c₁＝c₂＝0

$c_3=\frac{1}{3!}(-\frac{{\mathrm{\&beta;}}_3}{{\mathrm{\&delta;}}_3})$

$c_3=\frac{1}{3!}(-\frac{{\mathrm{\&beta;}}_3}{{\mathrm{\&delta;}}_3})$

$c_4=\frac{1}{4!}(-\frac{{\mathrm{\&beta;}}_4}{{\mathrm{\&delta;}}_4})$

$c_5=\frac{1}{5!}(-\frac{{\mathrm{\&beta;}}_5}{{\mathrm{\&delta;}}_5})$

$c_6=\frac{1}{6!}(-\frac{{\mathrm{\&beta;}}_6}{{\mathrm{\&delta;}}_6})$

将电压幅值和相角的各阶中心距代入，即可求解出电压幅值和相角的分布函数。

下面结合附图具体说明：

所述步骤1提供系统的各项分布数据，具体方法为：

假设风速服从三参数的威布尔分布，如图1所示，则风速的分布函数为：

$F_{\mathrm{Weibull}}\left(v\right)=1-\exp [-{\left(\frac{v-v_0}{c}\right)}^k]$

式中，v为风速，k、c、v₀分别代表威布尔分布的三个参数。其中k代表形状参数，它反映了风速的分布特点，c代表尺度参数，反映的是该地区的平均风速大小，v₀代表位置参数。

根据给定风机的切入、切出及额定风速，计算风机的输出功率，变化关系如图2所示，表达式为：

$P_W=\left\{\begin{array}{cc}0,& \\ k_{1}v+k_2,& v_{\mathrm{ci}}<v{\&le;v}_r\\ P_r,& v_r<v\&le;v_{\mathrm{co}}\\ 0,& v>v_{\mathrm{co}}\end{array}\right.$

式中：k₁＝P_r/(v_r-v_ci)；k₂＝-k₁v_ci；P_r为风力发电机的额定功率；v_ci为切入风速；v_r为额定风速；v_co为切出风速。

上述数据一般由风电场直接给出。

根据负荷功率的历史数据，计算期望和方差，通用公式如下：

如负荷功率数据x₁…x_n，对于服从正态分布的负荷，各数值的分布概率为p₁…p_n，计算负荷的期望m_x和方差δ_x：

$m_x=\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i\right)/n$

${\mathrm{\&delta;}}_x=\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_i-m_x)}^2\right)/n$

同理，将负荷风电场输出功率数据P_W代入，求解相应的期望与方差。

所述步骤2为基于RX模型的含风电场系统潮流计算，具体方法为：

图3中是考虑滑差的异步机RX模型。其中r₁+jx₁为定子阻抗，r₂+jx₂为转子阻抗，r_m+jx_m为激磁阻抗。自然风吹动风机叶片，将风能转化为机械能，由此获得的机械功率扣除掉机械损耗和杂散损耗后即为传递到异步发电机转子上的机械功率P_Ω，在等效电路中对应转子回路可变电阻r₂(1-s)/s上的功率。扣除转子铜耗P_Cu2和铁芯损耗P_Fe,得到输入定子绕组的电磁功率P_m，再扣除定子铜耗P_Cu1即得到注入电网的电功率P_e。图3中，由于x_m＞＞x₁，忽略定子电阻r₁及铁芯损耗P_Fe，可将励磁支路移至电路首端，得到异步发电机Γ型等值电路，如图4所示。

由图4所示的电路关系可得注入电网的电功率为：

$p_e=\frac{{-U}^2{r}_2/s}{{(r_2/s)}^2+x_k^2}$

式中，x_k＝x₁+x₂。

异步机吸收的无功与发出的有功之间的关系为：

$Q_e=\frac{r_2^2+x_k{(x_k+x_m)s}^2}{r_2{x}_{m}s}{P}_e$

上述方程构成了考虑滑差的异步机RX模型。

因此，在上述模型基础上，当节点i连接有风电场，与节点i对应的潮流方程为：

$\left\{\begin{array}{c}{P}_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)-P_{\mathrm{Li}}-V_i\underset{j\&Element;i}{\mathrm{\&Sigma;}}{V}_j(G_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})=0\\ Q_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)-Q_{\mathrm{Li}}-V_i\underset{j\&Element;i}{\mathrm{\&Sigma;}}{V}_j(G_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}})=0\\ P_{\mathrm{mi}}-P_{\mathrm{ei}}(V_i,{\mathrm{\&theta;}}_i,s_i)=0\end{array}\right.$

式中，P_ei、Q_ei分别表示风机有功功率和无功功率；P_Li、Q_Li分别表示负荷有功功率和无功功率；V_i、θ_i为节点i的电压幅值和相角；G_ij、B_ij分别为节点i、j间的电导与电纳值；θ_ij为节点i、j间的相角差；s_i为异步机滑差。

在上述潮流方程中，将未知量V_i、θ_i移至等式右侧，进行方程求解。

即，在已知风机有功功率和无功功率P_ei、Q_ei，负荷有功功率和无功功率P_Li、Q_Li，异步机滑差s_i，节点导纳G_ij、B_ij的情况下，根据上述潮流方程，计算风电场节点电压幅值和相角V_i、θ_i，通用表达式如下：

W＝f(X)

式中，W为节点注入向量，包括节点注入有功功率及无功功率P_ei、Q_ei，为已知。X为节点状态变量，包括节点的电压幅值和相角V_i、θ_i，为未知。

在概率潮流计算中，节点注入量是随机变量，因此可以将其表示为：

W＝W₀+ΔW

式中，W₀为节点注入向量W的期望值，ΔW为节点注入向量W的随机扰动。

同理可将状态变量写成

X＝X₀+ΔX

式中，X₀是状态变量X的期望值，ΔX为状态变量X的的随机变量。

进行泰勒级数展开，忽略高次项，得到：

W＝W₀+ΔW＝f(X₀+ΔX)＝f(X₀)+J₀ΔX

其中：

W₀＝f(X₀)

可以通过牛顿法潮流计算解出系统的状态变量X₀，J₀为牛顿法潮流计算最后一次迭代使用的雅克比矩阵。

ΔW＝J₀ΔX

因此有：

$\mathrm{\&Delta;X}={J_0}^{-1}\mathrm{\&Delta;W}=S_0\mathrm{\&Delta;W}$

式中，S₀为J₀的逆矩阵，称为灵敏度矩阵。

通过上述变换，计算得到灵敏度矩阵S₀各元素。

潮流计算及灵敏度矩阵的详细表述方法见《现代电力系统分析》（2003年，科学出版社）。

所述步骤3为各计算节点注入量的各阶半不变量，统一计算方法具体为：

如负荷功率数据x₁…x_n，计算其k阶中心矩，并记作β_k，计算方法如下：

${\mathrm{\&beta;}}_k=\underset{i=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(x_i-m_x)}^k{p}_i$

各阶半不变量可以由各自的各阶中心矩算得，通用计算方式如下：

γ₁＝m_x

${\mathrm{\&gamma;}}_2={{\mathrm{\&delta;}}_x}^2$

γ₃＝β₃

${\mathrm{\&gamma;}}_4={\mathrm{\&beta;}}_4-3{\mathrm{\&beta;}}_2^2$

γ₅＝β₅-10β₂β₃

${\mathrm{\&gamma;}}_6={\mathrm{\&beta;}}_6-15{\mathrm{\&beta;}}_2{\mathrm{\&beta;}}_4-10{\mathrm{\&beta;}}_3^2+30{\mathrm{\&beta;}}_2^3$

式中，γ_i代表i阶半不变量。

根据上述统一公式，由各变量的已知数据，计算节点i注入有功功率和无功功率的n阶半不变量

风电场节点的注入有功功率和无功功率的n阶半不变量

滑差及机械功率与电磁功率不平衡量的n阶半不变量

所述步骤4为由节点注入功率的各阶半不变量分别求取状态变量的各阶半不变量，具体步骤为：

根据下式，由注入量的各阶半不变量求得状态量

的半不变量：

$\left[\begin{array}{c}{\mathrm{\&Delta;\&theta;}}_1^{\left(n\right)}\\ {\mathrm{\&Delta;v}}_1^{\left(n\right)}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;\&theta;}}_n^{\left(n\right)}\\ {\mathrm{\&Delta;v}}_n^{\left(n\right)}\\ {\mathrm{\&Delta;s}}_n^{\left(n\right)}\end{array}\right]=\left[\begin{array}{ccccc}{S}_{11}^n& S_{12}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{\mathrm{1,2}n}^n& S_{\mathrm{1,2}n+1}^n\\ S_{21}^n& S_{22}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{\mathrm{2,2}n}^n& S_{\mathrm{2,2}n+1}^n\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;\\ S_{2n-\mathrm{1,1}}^n& S_{2n-\mathrm{1,2}}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n-\mathrm{1,2}n}^n& S_{2n-\mathrm{1,2}n+1}^n\\ S_{2n,1}^n& S_{2n,2}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n,2n}^n& S_{2n,2n+1}^n\\ S_{2n+\mathrm{1,1}}^n& S_{2n+1,2}^n& \&CenterDot;\&CenterDot;\&CenterDot;& S_{2n+\mathrm{1,2}n}^n& S_{2n+\mathrm{1,2}n+1}^n\end{array}\right]\left[\begin{array}{c}{\mathrm{\&Delta;P}}_1^{\left(n\right)}\\ {\mathrm{\&Delta;Q}}_1^{\left(n\right)}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}\\ {\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}\\ {\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}\end{array}\right]$

式中，

表示节点i的电压幅值和相角的n阶半不变量，

表示滑差及机械功率与电磁功率不平衡量的n阶半不变量；

表示灵敏度矩阵第i行第j列元素的n次方。

通过上式推导可得，状态量的半不变量可通过下式代数表达式计算得到：

${\mathrm{\&Delta;\&theta;}}_i^{\left(n\right)}=S_{2i-\mathrm{1,1}}^n{\mathrm{\&Delta;P}}_1^{\left(n\right)}+S_{2i-\mathrm{1,2}}^n{\mathrm{\&Delta;Q}}_1^{\left(n\right)}+\&CenterDot;\&CenterDot;\&CenterDot;+\&CenterDot;\&CenterDot;\&CenterDot;S_{2i-\mathrm{1,2}n-1}^n{\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}+S_{2i-\mathrm{1,2}n}^n{\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}+S_{2i-\mathrm{1,2}n+1}^n{\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}$

${\mathrm{\&Delta;v}}_i^{\left(n\right)}=S_{2i,1}^n{\mathrm{\&Delta;P}}_1^{\left(n\right)}+S_{2i,2}^n{\mathrm{\&Delta;Q}}_1^{\left(n\right)}+\&CenterDot;\&CenterDot;\&CenterDot;+\&CenterDot;\&CenterDot;\&CenterDot;S_{2i,2n-1}^n{\mathrm{\&Delta;P}}_{\mathrm{wind}}^{\left(n\right)}+S_{2i,2n}^n{\mathrm{\&Delta;Q}}_{\mathrm{wind}}^{\left(n\right)}+S_{2i,2n+1}^n{\mathrm{\&Delta;P}}_{\mathrm{M\&Omega;}}^{\left(n\right)}$

所述步骤5根据Gram-Charlier级数展开求得状态量的分布函数。具体步骤为：

通过节点i的电压幅值和相角的n阶半不变量

反推其各阶中心距，其统一表达式为：

γ₁＝m_x

${\mathrm{\&gamma;}}_2={{\mathrm{\&delta;}}_x}^2$

γ₃＝β₃

${\mathrm{\&gamma;}}_4={\mathrm{\&beta;}}_4-3{\mathrm{\&beta;}}_2^2$

γ₅＝β₅-10β₂β₃

${\mathrm{\&gamma;}}_6={\mathrm{\&beta;}}_6-15{\mathrm{\&beta;}}_2{\mathrm{\&beta;}}_4-10{\mathrm{\&beta;}}_3^2+30{\mathrm{\&beta;}}_2^3$

分布函数的统一表达式为:

F(ξ)＝c₀Φ(ξ)+c₁Φ′(ξ)+c₂Φ′′(ξ)+c₃Φ⁽³⁾(ξ)+c₄Φ⁽⁴⁾(ξ)+c₅Φ⁽⁵⁾(ξ)+c₆Φ⁽⁶⁾(ξ)式中，Φ(ξ)为标准正态分布函数表达式，各系数由下式求得

c₀＝1

c₁＝c₂＝0

$c_3=\frac{1}{3!}(-\frac{{\mathrm{\&beta;}}_3}{{\mathrm{\&delta;}}_3})$

$c_3=\frac{1}{3!}(-\frac{{\mathrm{\&beta;}}_3}{{\mathrm{\&delta;}}_3})$

$c_4=\frac{1}{4!}(-\frac{{\mathrm{\&beta;}}_4}{{\mathrm{\&delta;}}_4})$

$c_5=\frac{1}{5!}(-\frac{{\mathrm{\&beta;}}_5}{{\mathrm{\&delta;}}_5})$

$c_6=\frac{1}{6!}(-\frac{{\mathrm{\&beta;}}_6}{{\mathrm{\&delta;}}_6})$

将电压幅值和相角的各阶中心距代入，即可求解出电压幅值和相角的分布函数。

结果验证：为了测试本发明所提方法的有效性，应用本发明方法对IEEE-14和IEEE-30节点系统进行了仿真验证。

分别在两标准系统的14号、29号节点接入风电场。风电场单机容量为600kW,单台异步风力发电机的额定功率因数为0.89。风电场参数为：风电场空气密度为ρ＝1.2245kg/m³，风力机的扫掠面积为A＝1840m²。风电机组的切入风速、切出风速、额定风速均相同，分别为4m/s、25m/s、15m/s。异步风力发电机定子阻抗为0.00661+j0.07923Ω，转子阻抗为0.00296+j0.08654Ω，激磁电抗为2.81461Ω。本发明考虑负荷和风速的随机性，假设负荷服从正态分布，风速服从三参数的威布尔分布，风速的威布尔分布的三参数如下：v₀＝3；k＝3.97；c＝10.7。设所有风电机组型号相同，将n台风电机的机群等值为一台风力机，等值机参数如下：

$r_{\left(n\right)}=\frac{r}{n};$ $x_{\left(n\right)}=\frac{x}{n};$ s_(n)＝s

式中，r和x分别代表等值前单台风机的电阻和电抗参数，包括定子侧、转子侧和励磁，下标带(n)的表示相应的等值机参数。

利用本发明方法对IEEE-14、30系统进行仿真计算，依次加入20台、30台、40台风力发电机，IEEE-14节点系统风电场接入点处14号节点电压的分布函数如图5-图7所示，IEEE-30节点系统风电场接入点处29号节点电压的分布函数如图8-图10所示。

通过对含风电场的IEEE-14、302个系统进行概率潮流计算，并和蒙特卡罗法进行对比，结果表明本发明所提方法是正确和有效的。

Claims (6)

1.一种含风电场电力系统概率潮流的计算方法，包括如下步骤： 

步骤1：计算风电场功率以及负荷功率的概率分布； 

步骤2：用牛顿法进行含风电场电力系统确定性潮流的计算，求出灵敏度矩阵S₀； 

步骤3：计算各节点注入向量的各阶半不变量； 

步骤4：根据注入量的各阶半不变量分别求解状态变量的各阶半不变量； 

步骤5：根据Gram-Charlier级数展开求解概率密度函数和累积分布函数。 

2.根据权利要求1所述一种含风电场电力系统概率潮流的计算方法，其特征在于，所述步骤1包括： 

假设风速服从三参数的威布尔分布，则风速的分布函数F_Weibull(v)为： 

式中，v为风速，k、c、v₀分别代表威布尔分布的三个参数，其中k代表形状参数，c代表尺度参数，v₀代表位置参数； 

根据给定风机的切入、切出及额定风速，计算风机的输出功率P_W如下式： 

式中：k₁＝P_r/(v_r-v_ci)；k₂＝-k₁v_ci，P_r为风力发电机的额定功率，v_ci为切入风速，v_r为额定风速，v_co为切出风速； 

分别根据负荷功率数据x₁、x₂、…、x_n和风机的输出功率P_W的历史数据，计算期望和方差。 

3.根据权利要求1所述一种含风电场电力系统概率潮流的计算方法，其特征在于，所述步骤2包括： 

假设节点i连接有风电场，则与节点i对应的潮流方程为： 

式中，P_ei(V_i,θ_i,s_i)、Q_ei(V_i,θ_i,s_i)分别表示与变量V_i、θ_i和s_i有关的风机有功功率和无功功率；P_Li、Q_Li分别表示负荷有功功率和无功功率；V_i、θ_i分别为节点i的电压幅值和相角；G_ij、Bi_j分别为节点i、j间的电导与电纳值；θ_ij为节点i、j间的相角差；s_i为异步机滑差；P_mi为风力发电机的机械功率。 

求解上述潮流方程：在已知风机有功功率和无功功率P_ei(V_i,θ_i,s_i)、Q_ei(V_i,θ_i,s_i)，负荷有功功率和无功功率P_Li、Q_Li，异步机滑差s_i，G_ij、B_ij的情况下，根据上述潮流方程，计算风电场节点电压幅值V_i和相角θ_i，通用表达式如下： 

W＝f(X) 

式中，W为节点注入向量，包括P_ei(V_i,θ_i,s_i)和Q_ei(V_i,θ_i,s_i)，为已知；X为节点状态变量，包括节点的电压幅值V_i和相角θ_i，为未知； 

在概率潮流计算中，节点注入向量是随机变量，因此将其表示为： 

W＝W₀+ΔW 

式中，W₀为节点注入向量W的期望值，ΔW为节点注入向量W的随机扰动； 

同理将状态变量写成 

X＝X₀+ΔX 

式中，X₀是状态变量X的期望值，ΔX为状态变量X的随机变量； 

进行泰勒级数展开，忽略高次项，得到： 

W＝W₀+ΔW＝f(X₀+ΔX)＝f(X₀)+J₀ΔX 

其中： 

W₀＝f(X₀) 

通过牛顿法潮流计算解出系统的状态变量X₀，J₀为牛顿法潮流计算最后一次迭代使用的雅克比矩阵； 

ΔW＝J₀ΔX 

因此有： 

ΔX＝J₀ ^-1ΔW＝S₀ΔW 

式中，S₀为J₀的逆矩阵，称为灵敏度矩阵。 

4.根据权利要求1所述一种含风电场电力系统概率潮流的计算方法，其特征在于，所述步骤3包括： 

分别计算负荷功率数据x₁、x₂、…、x_n的k阶中心矩，并记作β_k，计算方法如下： 

式中，m_x为负荷功率的期望，p_i为x₁、x₂、…、x_n的分布概率，其中i=1、2、……、n，各阶半不变量由各自的各阶中心矩算得，通用计算方式如下： 

γ₁＝m_x

γ₂＝δ_x ²

γ₃＝β₃

γ₅＝β₅-10β₂β₃

式中，γ_i代表i阶半不变量，δ_x为负荷功率的方差； 

根据上述公式，由各变量的已知数据，分别计算节点i注入有功功率和无功功率的n阶半不变量

和

风电场节点的注入有功功率和无功功率的n阶半不变量

和

分别计算滑差与电磁功率不平衡量的n阶半不变量和机械功率与电磁功率不平衡量的n阶半不变量。 

5.根据权利要求1所述一种含风电场电力系统概率潮流的计算方法，其特征在于，所述步骤4包括： 

根据下式，由注入量的各阶半不变量求得状态量

的半不变量： 

式中，

和

表示节点i的电压幅值和相角的n阶半不变量，

和

分别表示滑差与电磁功率不平衡量的n阶半不变量和机械功率与电磁功率不平衡量的n阶半不变量；

表示灵敏度矩阵第i行第j列元素的n次方； 

通过上式推导，状态变量的半不变量通过下式计算得到： 

6.根据权利要求4所述一种含风电场电力系统概率潮流的计算方法，其特征在于，所述步骤5包括： 

通过节点i的电压幅值和相角的n阶半不变量

和

反推其各阶中心距，其统一表达式为： 

γ₁＝m_x

γ₂＝δ_x ²

γ₃＝β₃

γ₅＝β₅-10β₂β₃

分布函数F(ξ)的统一表达式为: 

F(ξ)＝c₀Φ(ξ)+c₁Φ′(ξ)+c₂Φ′′(ξ)+c₃Φ⁽³⁾(ξ)+c₄Φ⁽⁴⁾(ξ)+c₅Φ⁽⁵⁾(ξ)+c₆Φ⁽⁶⁾(ξ)式中，Φ(ξ)为标准正态分布函数表达式，各系数由下式求得 

c₀＝1 

c₁＝c₂＝0 

将电压幅值和相角的各阶中心距代入，即求解出电压幅值和相角的分布函数。 

Images