CN105207204A - 一种计及一次调频不确定性的概率潮流分析方法 - Google Patents

Abstract

本发明公开提出了一种计及一次调频不确定性的概率潮流分析方法，构造计及一次调频特性的动态潮流模型，应用点估计方法处理模型中的不确定问题，采用牛拉法求解确定性潮流问题，从而求出扰动不确定情况下，计及系统一次调频特性及其系数不确定性时的系统频率、节点电压等状态变量的统计特征值。在计及一次调频特性的动态潮流模型设计中，考虑发电机和负荷单位调节功率必然存在的不确定性，因此可以得到更接近实际情况的状态变量，为相关的分析计算提供更准确的数据。该方法可以更为有效准确地描述系统的运行状态，为电力系统进一步的分析计算提供基础。同时，点估计法具有速度快，低阶准确率高等优点，为本方法的实现提供了有效途径。

Description

一种计及一次调频不确定性的概率潮流分析方法

技术领域

本发明涉及电力系统稳态潮流计算技术领域，具体涉及一种计及一次调频不确定性的概率潮流分析方法。

背景技术

近年来随着化石能源的日益枯竭、环境污染问题的愈加严重以及新能源发电技术的不断发展，以风力发电为代表的具有清洁环保、可再生特点的新能源接入电网的规模逐渐扩大，其发电量占系统总发电量的比例也不断增加。但于此同时，这些新能源所固有的随机性和间歇性等特点也使得电网运行时的不确定性和波动性愈加明显。因此，如何有效地在潮流计算中考虑这些不确定性和扰动性，以获得更为准确的系统状态，显得尤为重要。潮流计算作一种获得系统节点电压、支路功率等状态信息的有效和常用方法，是计算最优潮流、状态估计等一系列电网分析的基础，具有重要的研究价值。

自1974年以来，概率潮流作为一种研究不确定条件下潮流分布的有效方法得到充分的发展和应用。概率潮流的基本原理是结合随机参数的特征值及其概率分布，通过多次潮流计算，从而获得系统状态参数的特征值及其概率分布。长期以来，概率潮流都建立在常规潮流模型的基础上，也就是认为当处于稳定状态的系统出现扰动时，所有不平衡功率全部由平衡机组承担，其余机组的出力和系统各负荷功率值均保持不变。然而在实际系统中，平衡机组容量有限，不可能满足系统任何情况下的功率波动。当系统中扰动过大而平衡机组又无法满足时，必将导致系统频率的波动，且由于系统固有的一次调频特性，非平衡机组的出力以及各负荷节点的功率值也将随之改变。此时，若仍采用基于常规潮流模型的概率潮流进行计算，则所得结果与实际运行状态差别较大，并且扰动越大，误差越大。因此，随着新能源大规模地接入电网，有必要在概率潮流中计及系统的一次调频特性，以得到更加接近实际的潮流分布。

系统的一次调频特性由发电机的静态频率特性和负荷的静态频率特性两部分组成。其中发电机是否投入一次调频以及频率偏差多大时进行频率调节不仅与机组容量有关，还与机组类型、死区以及调差参数R的设置等相关。但由于各种原因，这些参数的准确值一般很难获得。曾对多次事故后发电机组的数据进行分析，发现每次事故中均存在一些机组未能按照预定的功频特性调整其出力的情况。这是由于机组的死区与调差参数R未按照要求设置造成。由此可见，发电机是否参与一次调频存在随机性，另一方面，负荷的单位调节功率与负荷的类型直接相关，不同的负荷组成，其一次调频性能存在很大的差异。因此，在一次调频计算时计及其调节系数的随机性对于准确地描述系统状态具有非常重要的意义。

发明内容

本发明的目的是针对现有概率潮流研究方法的不足，提供了一种计及一次调频不确定性的概率潮流分析方法，并采用点估计法进行求解。该方法可以更有效准确的描述系统的运行状态，为电力系统进一步的分析计算提供基础。同时，点估计法具有速度快，低阶准确率高等优点，为本方法的实现提供了有效途径。

为实现本发明目的而采用的技术方案是这样的，构造计及一次调频特性的动态潮流模型，应用点估计方法处理模型中的不确定问题，采用牛拉法求解确定性潮流问题，从而求出扰动不确定情况下，计及系统一次调频特性及其系数不确定性时的系统频率、节点电压等状态变量的统计特征值。

一种计及一次调频不确定性的概率潮流分析方法，包括以下具体步骤：

1)输入确定参数及随机变量的分布信息；

输入常规潮流计算所需要的基础信息、随机变量的分布信息和随机变量的数量。

所述输入常规潮流计算所需要的基础信息，包括：网络拓扑数据，支路阻抗，母线并联电导电纳，节点类型，节点负荷功率，PV节点上发电机出力数据，PQ节点上发电机出力数据，发电机的额定容量及功率因素，发电机的调差系数标幺值R_g和系统负荷扰动功率预测值。

所述随机变量的分布信息，其中随机变量包括：节点负荷扰动预测误差ΔP_{L_pre}、负荷的单位调节功率k_L、发电机的单位调节功率k_G。所述预测误差ΔP_{L_pre}服从正态分布，确定其正态分布的期望和方差；所述发电机的单位调节功率k_G服从二项分布，确定发电机单位调节功率标幺值的两点取值(即K_Gg＝1/R_g和K_Gg＝0)及其对应的概率；所述负荷的单位调节功率k_L服从均匀分布，确定其标幺值k_Lg的分布区间[a',b']。

所述随机变量的数量，为输入系统中随机变量的数量，总个数m，其中随机变量ΔP_{L_pre}的个数为m_p,随机变量k_L的个数为m_kL，随机变量k_G的个数为m_kG。

2)计算基准状态下的潮流分布；

结合步骤1获得的所述常规潮流基础信息，建立如下常规潮流模型：

$\left\{\begin{array}{c}{P}_{Gi}-P_{Li}^0=U_i\underset{j=1}{\overset{j=n}{\Σ}}{U}_j(G_{ij}{\mathrm{cos\δ}}_{ij}+B_{ij}{\mathrm{sin\δ}}_{ij})\\ Q_{Gi}-Q_{Li}^0=U_i\underset{j=1}{\overset{j=n}{\Σ}}{U}_j(G_{ij}{\mathrm{sin\δ}}_{ij}-B_{ij}{\mathrm{cos\δ}}_{ij})\end{array}\right.---\left(1\right)$

公式1是基于极坐标的常规潮流方程，其中P_Gi,Q_Gi分别表示节点i所接发电机的有功出力和无功出力，若该节点无发电机，则P_Gi和Q_Gi均取0。分别表示基准状态下，节点i上负荷吸收的有功功率和无功功率，若该节点无负荷，则和均取0。U_i,U_j分别表示节点i和节点j的电压幅值，δ_ij＝θ_i-θ_j表示节点i、j间电压的相角差，其中θ_i和θ_j分别表示节点i和节点j电压的相角。G_ij,B_ij分别表示节点i、j间支路的电导和电纳，若节点i、j间不存在支路则对应值为零。

运用牛拉法求解上述方程，得到基准状态下发电机有功出力无功出力节点电压幅值相角

3)构建计及系统一次调频特性的动态潮流模型；

结合步骤1获得的数据以及步骤2得到的基准状态变量，构建负荷扰动时，计及负荷和发电机一次调频作用的动态潮流模型如下：

$\left\{\begin{array}{c}{P}_{Gi}-P_{Li}=U_i\underset{j=1}{\overset{j=n}{\Σ}}{U}_j(G_{ij}{\mathrm{cos\δ}}_{ij}+B_{ij}{\mathrm{sin\δ}}_{ij})\\ Q_{Gi}^0-Q_{L1}^0=U_i\underset{j=1}{\overset{j=n}{\Σ}}{U}_j(G_{ij}{\mathrm{sin\δ}}_{ij}-B_{ij}{\mathrm{cos\δ}}_{ij})\\ P_{Gi}=P_{Gi}^0-k_{Gi}\mathrm{\Δ}f\\ P_{Li}=P_{Li}^{\′}+k_{Li}\mathrm{\Δ}f=(P_{Li}^0+{\mathrm{\ΔP}}_{Li})+k_{Li}\mathrm{\Δ}f\end{array}\right.---\left(2\right)$

公式2中，P_Gi,P_Li分别表示系统发生扰动后，计及一次调频作用时节点i上发电机的实际有功出力和负荷吸收的实际有功功率。k_Gi,k_Li分别表示节点i上发电机和负荷的单位调节功率，Δf＝f-f₀为系统的频率偏移，其中f₀和f分别表示基准状态的系统频率(50Hz)和扰动后的系统频率。P_L'_i表示扰动发生瞬间节点i上的有功，ΔP_Li分别表示扰动前基准状态节点i上的有功负荷及实际产生的功率扰动，分别表示步骤2求得的基准状态下节点i上发电机的有功出力和无功出力。其余变量的含义与步骤2相同。

由于本发明仅关注系统单位调节功率不确定性对频率、发电机有功出力等状态变量的影响，且考虑到频率与有功强耦合，与无功弱耦合，故本发明中不计及无功调整。

同时计及系统的一次调频特性，系统中的不平衡功率不再仅由平衡节点消纳，而是由系统中所有发电机和有功负荷一起分担，因此不再有平衡节点的概念，且认为原平衡节点的无功出力仍能保证该点电压幅值恒定不变。因此，相对常规潮流模型，动态潮流模型需要增加频率变量以及原平衡节点的有功不平衡方程。

4)根据分布信息，计算随机变量所需的统计特征值；

对概率密度为f_X(x)的变量X,求取下列统计特征值：

期望μ_X：

${\mathrm{\μ}}_X=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{\mathrm{xf}}_X\left(x\right)dx---\left(3\right)$

标准差σ_X：

${\mathrm{\σ}}_X=\sqrt{\left(\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{(x-{\mathrm{\μ}}_X)}^2{f}_X\right(x\left)dx\right)}---\left(4\right)$

X的j阶中心距M_j(X)：

$M_j\left(X\right)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{(x-{\mathrm{\μ}}_X)}^j{f}_X\left(x\right)dx---\left(5\right)$

定义λ_{X_j}表示随机变量X的j阶中心距与标准差σ_X的j次方的比值：

${\mathrm{\λ}}_{X\_j}=\frac{M_j\left(X\right)}{{\mathrm{\σ}}_X^j}---\left(6\right)$

其中，λ_{X_1}＝0，λ_{X_2}＝1，λ_{X_3}和λ_{X_4}分别表示随机变量X的偏度系数和峰度系数。

结合公式3-6，并根据步骤1得到的随机变量分布信息，计算各随机变量的期望、均方差、偏度系数和峰度系数。包括以下过程：

I)负荷扰动预测误差；

实际的扰动功率由两部分组成：

ΔP_L＝P_{L_pre}+ΔP_{L_pre}(7)

公式7中，ΔP_L、P_{L_pre}分别表示实际的负荷扰动及其预测值，ΔP_{L_pre}表示预测误差，且认为其服从标准差为σ的零期望正态分布。

用变量X1表示预测误差ΔP_{L_pre}，其各统计特征值求取如下：

计算扰动预测误差ΔP_{L_pre}的偏度系数λ_{X1_3}：

${\mathrm{\λ}}_{X1\_3}=\frac{M_3\left(X1\right)}{{\mathrm{\σ}}_{X1}^3}=\frac{0}{{\mathrm{\σ}}_{X1}^3}=0---\left(8\right)$

计算扰动预测误差ΔP_{L_pre}的峰度系数λ_{X1_4}：

${\mathrm{\λ}}_{X1\_4}=\frac{M_4\left(X1\right)}{{\mathrm{\σ}}_{X1}^4}=\frac{3{\mathrm{\σ}}_{X1}^4}{{\mathrm{\σ}}_{X1}^4}=3---\left(9\right)$

II)发电机的单位调节功率；

机组的单位调节功率k_G服从二项分布，参与一次调频的概率为p，对应的单位调节功率的值为1/R,其中R是机组的调差系数，不参与一次调频的概率为q(q＝1-p)，对应的单位调节功率的值为0。

由于步骤1中输入的是调差系数的标幺值，因此需要先将其转换为有名值，公式如下：

$R=R_g\frac{f_N}{P_{GN}}---\left(10\right)$

其中f_N＝50Hz，为系统额定频率。P_GN是机组的额定有功功率，通过公式11求取：

P_GN＝S_GNcosα(11)

其中S_GN机组的额定容量，α为功率因素。

用变量X2表示发电机的单位调节功率k_G，其各统计特征值求取如下：

计算发电机的单位调节效率k_G的期望μ_X2：

${\mathrm{\μ}}_{X2}=p\frac{1}{R}+q\×0=\frac{p}{R}---\left(12\right)$

计算发电机的单位调节效率k_G的均方差σ_X2：

${\mathrm{\σ}}_{X2}=\sqrt{p{(\frac{1}{R}-\frac{p}{R})}^2+q{(0-\frac{p}{R})}^2}=\frac{\sqrt{pq}}{R}---\left(13\right)$

计算发电机的单位调节效率k_G的偏度系数λ_{X2_3}：

${\mathrm{\λ}}_{X2\_3}=\frac{M_3\left(X2\right)}{{\mathrm{\σ}}_{X2}^3}=\frac{p{(1/R-p/R)}^3+q{(0-p/R)}^3}{{(\sqrt{pq}/R)}^3}=\frac{q^2-p^2}{\sqrt{pq}}---\left(14\right)$

计算发电机的单位调节效率k_G的峰度系数λ_{X3_4}：

${\mathrm{\λ}}_{X3\_4}=\frac{M_4\left(X3\right)}{{\mathrm{\σ}}_{X3}^4}=\frac{p{(1/R-p/R)}^4+q{(0-p/R)}^4}{{(\sqrt{pq}/R)}^4}=\frac{q^3+p^3}{pq}---\left(15\right)$

III)负荷的单位调节功率；

负荷的单位调节率标幺值k_Lg服从均匀分布，区间为[a',b']。与步骤II一样，需要先将标幺值区间[a',b']转换为有名值区间[a,b]，按照公式16：

$k_L=k_{Lg}\frac{P_{LN}}{f_N}---\left(16\right)$

其中P_LN表示负荷的额定有功，认为其等于基准状态下负荷实际吸收的有功，即步骤1中输入的负荷功率；f_N＝50Hz，为系统额定频率。

用变量X3表示负荷的单位调节功率k_L，其各统计特征值求取如下：

利用公式17计算负荷的单位调节功率k_L的期望μ_X3：

${\mathrm{\μ}}_{X3}=\frac{b-a}{2}---\left(17\right)$

利用公式18计算负荷的单位调节功率k_L的均方差σ_X3：

${\mathrm{\σ}}_{X3}=\sqrt{\underset{a}{\overset{b}{\∫}}{(x-{\mathrm{\μ}}_{X3})}^2\frac{1}{b-a}dx}=\frac{b-a}{2\sqrt{3}}---\left(18\right)$

利用公式19计算负荷的单位调节功率k_L的偏度系数λ_{X3_3}：

${\mathrm{\λ}}_{X3\_3}=\frac{M_3\left(X3\right)}{{\mathrm{\σ}}_{X3}^3}=\frac{\underset{a}{\overset{b}{\∫}}{(x-{\mathrm{\μ}}_{X3})}^3\frac{1}{b-a}dx}{{\mathrm{\σ}}_{X3}^3}=0---\left(19\right)$

利用公式20计算负荷的单位调节功率k_L的峰度系数λ_{X3_4}：

${\mathrm{\λ}}_{X3\_4}=\frac{M_4\left(X3\right)}{{\mathrm{\σ}}_{X3}^4}=\frac{\underset{a}{\overset{b}{\∫}}{(x-{\mathrm{\μ}}_{X3})}^4\frac{1}{b-a}dx}{{\mathrm{\σ}}_{X3}^4}=\frac{9}{5}---\left(20\right)$

5)根据2m+1点估计法构造随机变量的样本点，确定样本值及其权重系数；

根据点估计原理，对系统中的所有随机变量z_i(z_i∈{z₁,z₂,…z_m})分别按照公式21构造三个样本值z_i,1,z_i,2,z_i,3：

z_i,j＝μ_i+ξ_i,jσ_ii＝1,…m；j＝1,2,3(21)

其中ξ_i,j表示第i个随机变量的第j个位置系数，通过下式求取：

$\left\{\begin{array}{c}{\mathrm{\ξ}}_{i,1}=\frac{{\mathrm{\λ}}_{i,3}}{2}+\sqrt{{\mathrm{\λ}}_{i,4}-\frac{3}{4}{\mathrm{\λ}}_{i,3}^2}\\ {\mathrm{\ξ}}_{i,2}=\frac{{\mathrm{\λ}}_{i,3}}{2}-\sqrt{{\mathrm{\λ}}_{i,4}-\frac{3}{4}{\mathrm{\λ}}_{i,3}^2}\\ {\mathrm{\ξ}}_{i,3}=0\end{array}\right.---\left(22\right)$

公式22中，λ_i3和λ_i4分别表示步骤4中求取的第i个随机变量的偏度系数和峰度系数。

基于公式22得到的参数，求取公式22构造的三个样本值对应的权重系数w_i,1,w_i,2,w_i,3，

$\left\{\begin{array}{c}{w}_{i,1}=\frac{1}{{\mathrm{\ξ}}_{i1}^2-{\mathrm{\ξ}}_{i1}{\mathrm{\ξ}}_{i2}}\\ w_{i,2}=\frac{1}{{\mathrm{\ξ}}_{i2}^2-{\mathrm{\ξ}}_{i1}{\mathrm{\ξ}}_{i2}}\\ w_{i,3}=\frac{1}{m}-\frac{1}{{\mathrm{\λ}}_{i,4}-{\mathrm{\λ}}_{i,3}^2}\end{array}\right.---\left(23\right)$

其中m表示系统中随机变量的个数。

然后，基于得到的负荷扰动预测误差的三个样本点，按照步骤4中的公式7，计算得到负荷扰动的样本。

完成以上步骤后，接着构造系统样本点。构造随机变量样本值z_i,j,(i＝1,…m；j＝1,2,3)对应的系统样本点时，随机变量z_i＝z_i,j，其余变量取均值，即样本点Z(i,j)＝(μ₁,…,z_i,j,…μ_m)，总共需要构造3m个系统样本点。但由公式22可知ξ_i3＝0(i＝1,…,m)，因此有Z(1,3)＝Z(2,3)＝…＝Z(m,3)＝(μ₁,…,μ_i,…μ_m)，故只需计算一个样本，减少m-1次计算，共计只需2m+1次确定潮流计算。

6)参数初始化；

完成步骤5后，给定潮流计算所需的初始参数。设置最大迭代次数Tmax，收敛精度ε＝10^-5～10^-3。根据步骤2得到的基准状态数据，给定发电机有功出力、无功出力初值分别为和节点电压幅值和相角初值分别为和系统频率初值f₀＝50Hz。

利用公式24计算节点i的注入无功：

$Q_{i\_input}^0=Q_{Gi}^0-Q_{Li}^0---\left(24\right)$

7)采用牛拉法计算各样本对应的系统状态变量值；

完成步骤6后，基于步骤5得到的由负荷扰动、发电机单位调节功率和负荷单位调节功率组成的系统样本和步骤3构造的动态潮流模型，分别计算3m个样本点对应的确定性潮流，获得该样本点对应的系统频率f，节点电压幅值U和相角θ，发电机的实际有功出力P_G和负荷吸收的实际有功功率P_L.。主要步骤包括：

i)设置迭代次数初值，计算初始不平衡量；

首先设置迭代次数初值t＝0。

然后计算基准状态发生扰动ΔP_Li瞬间的初始不平衡量和

$\left\{\begin{array}{cc}{\mathrm{\ΔP}}_i^0=P_{i\_input}^0-U_i^0\underset{j=1}{\overset{j=n}{\Σ}}{U}_j^0\left(G_{ij}{\mathrm{cos\δ}}_{ij}^0+B_{ij}{\mathrm{sin\δ}}_{ij}^0\right)& \left(i\∈\[pq;pv;ref\]\right)\\ {\mathrm{\ΔQ}}_i^0=Q_{i\_input}^0-U_i^0\underset{j=1}{\overset{j=n}{\Σ}}{U}_j^0\left(G_{ij}{\mathrm{sin\δ}}_{ij}^0-B_{ij}{\mathrm{cos\δ}}_{ij}^0\right)& \left(i\∈pq\right)\end{array}\right.---\left(25\right)$

公式25中为基准状态节点i电压相角和节点j电压相角之差，即：

${\mathrm{\δ}}_{ij}^0={\mathrm{\θ}}_i^0-{\mathrm{\θ}}_j^0---\left(26\right)$

表示扰动后节点i的注入有功，根据下式求取：

$\left\{\begin{array}{c}{P}_{i\_input}^0=P_{Gi}^0-P_{Li}^{\′0}\\ P_{Li}^{\′0}=P_{Li}^0+{\mathrm{\ΔP}}_{Li}=P_{Li}^0+P_{Li\_pre}+{\mathrm{\ΔP}}_{Li\_pre}\end{array}\right.---\left(27\right)$

其中和分别表示基准状态下节点i上发电机的有功出力和有功负荷；表示扰动瞬间节点i的有功负荷；P_{L_pre}和ΔP_{L_pre}分别表示扰动预测值和扰动预测误差。

此时系统的不平衡量为：

ii)计算雅克比矩阵；

基于步骤3提出的动态潮流模型，计算雅克比矩阵。

$J=\left[\begin{array}{ccc}H& N& M\\ K& L& T\\ W& O& V\end{array}\right]---\left(29\right)$

其中H为n_pq+n_pv阶的方阵：

$H_{ij}=\frac{\∂{\mathrm{\ΔP}}_i}{\∂{\mathrm{\θ}}_j}=\left\{\begin{array}{cc}{U}_i^t{U}_j^t(G_{ij}{\mathrm{sin\δ}}_{ij}^t-B_{ij}{\mathrm{cos\δ}}_{ij}^t)& (i\≠j)\\ -U_i^t\underset{k\≠i}{\underset{k=1}{\overset{n}{\Σ}}}{U}_k^t(G_{ik}{\mathrm{sin\δ}}_{ik}^t-B_{ik}{\mathrm{cos\δ}}_{ik}^t)& (i=j)\end{array}\right.,(i\∈\[pq;pv\];j\∈\[pq;pv\])---\left(30\right)$

N为(n_pq+n_pv)×n_pq阶的矩阵：

$N_{ij}=\frac{\∂{\mathrm{\ΔP}}_i}{\∂U_j}{U}_j=\left\{\begin{array}{cc}{U}_i^t{U}_j^t(G_{ij}{\mathrm{cos\δ}}_{ij}^t+B_{ij}{\mathrm{sin\δ}}_{ij}^t)& (i\≠j)\\ U_i^t\underset{k=1}{\overset{n}{\Σ}}{U}_k^t(G_{ik}{\mathrm{cos\δ}}_{ik}^t+B_{ik}{\mathrm{sin\δ}}_{ik}^t)+{\left(U_i^t\right)}^2{G}_{ii}& (i=j)\end{array}\right.,(i\∈\[pq;pv\];j\∈pq)---\left(31\right)$

M为(n_pq+n_pv)×1阶的矩阵：

$M_i=\frac{\∂{\mathrm{\ΔP}}_i}{\∂f}=(k_{Gi}+k_{Li}),(i\∈\[pq:pv\])---\left(32\right)$

公式32中若节点i上没有发电机则k_Gi＝0，若没有负荷则k_Li＝0；

K为n_pq×(n_pq+n_pv)阶的矩阵：

$K_{ij}=\frac{\∂{\mathrm{\ΔQ}}_i}{\∂{\mathrm{\θ}}_j}=\left\{\begin{array}{cc}-U_i^t{U}_j^t(G_{ij}{\mathrm{cos\δ}}_{ij}^t+B_{ij}{\mathrm{sin\δ}}_{ij}^t)& (i\≠j)\\ U_i^t\underset{k\≠i}{\underset{k=1}{\overset{n}{\Σ}}}{U}_k^t(G_{ik}{\mathrm{cos\δ}}_{ik}^t+B_{ik}{\mathrm{sin\δ}}_{ik}^t)& (i=j)\end{array}\right.,(i\∈pq;j\∈\[pq;pv\])---\left(33\right)$

L为n_pq×n_pq阶的矩阵：

$L_{ij}=\frac{\∂{\mathrm{\ΔQ}}_i}{\∂U_j}{U}_j=\{\begin{array}{cc}{U}_i^t{U}_j^t\left(G_{ij}{\mathrm{sin\δ}}_{ij}^t-B_{ij}{\mathrm{cos\δ}}_{ij}^t\right)& \left(i\≠j\right)\\ U_i^t\underset{k=1}{\overset{n}{\Σ}}{U}_k^t\left(G_{ik}{\mathrm{sin\δ}}_{ik}^t-B_{ik}{\mathrm{cos\δ}}_{ik}^t\right)-{\left(U_i^t\right)}^2{B}_{ii}& \left(i=j\right)\end{array},\left(i\∈pq;j\∈pq\right)---\left(34\right)$

T为n_pq×1阶的矩阵：

$T_i=\frac{\∂Q_i}{\∂f}=0,(i\∈pq)---\left(35\right)$

W为1×(n_pq+n_pv)阶的矩阵：

$W_j=\frac{\∂{\mathrm{\ΔP}}_i}{\∂{\mathrm{\θ}}_j}=U_i^t{U}_j^t(G_{ij}{\mathrm{sin\δ}}_{ij}^t-B_{ij}{\mathrm{cos\δ}}_{ij}^t),(i=ref,j\∈\[pq;pv\])---\left(36\right)$

O为1×n_pq阶的矩阵：

$O_j=\frac{\∂{\mathrm{\ΔP}}_i}{\∂U_j}{U}_j=U_i^t{U}_j^t(G_{ij}{\mathrm{cos\δ}}_{ij}^t+B_{ij}{\mathrm{sin\δ}}_{ij}^t),(i=ref;j\∈pq)---\left(37\right)$

V为1×1阶的矩阵：

$V=\frac{\∂{\mathrm{\ΔP}}_{ref}}{f}=k_{Gi}+k_{Li},(i=ref)---\left(38\right)$

其pq,pv分别表示PQ节点序号矩阵，PV节点序号矩阵，ref为平衡节点的序号。n,n_pv,n_pq,n_ref(n_ref＝1)分别表示系统总的节点个数，PQ节点个数，PV节点个数和平衡节点个数。

迭代次数自增，即t＝t+1。

iii)解修正方程，计算修正量；

基于步骤i和步骤ii，求解修正方程39，得到修正量：

ΔX＝-J^-1ΔF(39)

ΔX为变量的修正量如公式40：

对电压修正量进行公式41的变换：

${\mathrm{\ΔU}}_i=\frac{{\mathrm{\ΔU}}_i}{U_i}{U}_i^{t-1},(i\∈pq)---\left(41\right)$

其中表示上次迭代得到的节点i的电压幅值，初次迭代时指给定的电压幅值初值。

进一步，结合公式40和41，计算得到的修正量修正原有的变量，得到新值：

X^t＝X^t-1+ΔX(42)

iv)计算不平衡量并进行收敛判断；

按照公式25计算新值下的不平衡量ΔF：

判断不平衡量ΔF是否满足收敛判据：

max{|ΔF|}≤ε(44)

若收敛判据满足则迭代结束转到下一步；若收敛不满足且迭代次数t＜Tmax则跳转至步骤ii，继续迭代直至满足收敛条件；若判据不满足且t＝Tmax，则表示迭代不收敛，跳出循环，转至下一步。

v)结果输出；

判断此时t是否满足t＜Tmax，若满足则此时得到的结果X即为样本Z(i,j)对应的状态变量值；否则表示潮流迭代不收敛。

8)基于所得结果及其对应权重系数，计算变量的统计特征值；

基于步骤7所得的3m个样本对应的潮流计算结果以及步骤5计算得到的各样本点对应的权重系数，根据公式45-47计算变量的期望和方差。

$E\[x^k\]=\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{3}{\Σ}}{w}_{i,j}{\left(X\right(i,j\left)\right)}^k---\left(45\right)$

其中

X(1,3)＝…＝X(m,3)(46)

$\left\{\begin{array}{c}{\mathrm{\μ}}_x=E\[x\]\\ {\mathrm{\σ}}_x=\sqrt{E\[x^2\]-{{\mathrm{\μ}}_x}^2}\end{array}\right.---\left(47\right)$

其中X(i,j)表示样本Z(i,j)计算得到的变量x的值；μ_x和σ_x分别表示变量x的期望和均方差。变量X包含节点电压相角，电压幅值和频率，

本发明的技术效果是毋庸置疑的，采用上述技术方案后，主要有以下效果：

1)与现有考虑扰动的概率分析方法相比，本方案计及系统的一次调频作用，并考虑了发电机和负荷单位调节功率必然存在的不确定性，因此可以得到更接近实际情况的状态变量，为相关的分析计算提供更准确的数据。

2)本方案采用点估计方法，具有速度快，低阶精度高的特点，可以快速有效地进行计及不确定因素的概率潮流计算。

本方案可以广泛应用于扰动发生之前预测系统的潮流分布的情况，对于正确判断系统状态，并基于潮流分布的进一步相关分析具有非常重要的意义。

附图说明

图1为本发明方法的整体流程图；

图2为本发明方法中步骤7采用的牛拉法流程图。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

以IEEE-9节点系统为例，一种计及一次调频不确定性的概率潮流分析方法的具体步骤如下：

1)输入确定参数及随机变量的分布信息；

输入常规潮流计算所需要的基础信息，包括：网络拓扑数据、支路阻抗、母线并联电导电纳、节点类型、节点负荷功率，PV、PQ节点上发电机出力数据。

输入节点5、7和9有功功率扰动的预测值(单位：MW)；

P_{L_pre}＝[91012.5]

输入上述节点预测误差ΔP_{L_pre}的期望μ和均方根σ：

μ＝[000]

σ＝[0.911.25]

输入负荷单位调节功率标幺值k_Lg的取值范围：[0，3]。

输入发电机单位调节功率标幺值k_Gg的分布信息：参与调频的概率p＝0.9，调差系数R的标幺值为0.04，对应的发电机单位调节功率的标幺值为25；不参与调频的概率q＝0.1，对应的发电机单位调节功率的标幺值为0。输入发电机的额定容量(单位：MVA)S_GN＝[247.5，192.0，128.0]和功率因素α＝[1.0，0.85，0.85]。

输入总随机变量的个数m＝9，其中负荷扰动预测误差的个数mp＝3；负荷单位调节功率的个数为m_kL＝3；发电机单位调节功率的个数为m_kG＝3。

2)计算基准状态下的潮流分布，获得基准状态信息；

步骤1完成后，牛拉法计算基准状态下的发电机有功出力无功出力节点电压幅值相角潮流模型采用技术方案中的公式1。

根据IEEE-9节点系统的网络结构和参数信息，按照技术方案中公式1，运用牛拉法计算得到基准状态下潮流分布：

得到发电机有功出力(单位：MW)

$P_G^0={\left[\begin{array}{ccc}71.9547& 163& 85\end{array}\right]}^T$

得到无功出力(单位MVar)：

$Q_G^0={\left[\begin{array}{ccc}24.0690& 14.4601& -3.6490\end{array}\right]}^T$

得到节点电压幅值U⁰(标幺值)：

U⁰＝[1110.98700.97551.00340.98560.99620.9576]^T

得到相角θ⁰(单位：度)如下：

θ⁰＝[09.66874.7711-2.4066-4.01731.92560.62153.7991-4.3499]^T

3)构建计及系统一次调频特性的动态潮流模型；

按照技术方案中的公式2，构造计及系统一次调频特性的动态潮流模型。

4)根据分布信息，计算随机变量所需的统计特征值；

完成步骤2后，基于步骤1输入的随机变量分布信息，计算各随机变量的期望、方差、偏度系数和峰度系数。

I)负荷扰动预测误差；

根据步骤1输入的随机变量的分布信息，得到负荷扰动预测误差的期望μ_X1和方差σ_X1

μ_X1＝[000]

σ_X1＝[0.911.25]

按照技术方案中的公式8，计算得到偏度系数λ_{X1_3}：

λ_{X1_3}＝[000]^T

按照技术方案中的公式9，计算得到峰度系数λ_{X1_4}：

λ_{X1_4}＝[333]^T

II)发电机的单位调节功率；

根据输入信息，基于技术方案中的公式11计算得到机组的额定有功P_GN(单位MW)

P_GN＝[247.5163.2108.8]

按照技术方案中的公式10将调差系数R的标幺值转换为有名值(单位Hz/MW)：

R＝[0.00080.01230.0184]^T

按照技术方案中的公式12，计算得到k_G的期望μ_X2(单位MW/Hz)：

μ_X2＝[111.37573.4448.96]^T

按照技术方案中的公式13，计算得到k_G的均方差σ_X2：

σ_X2＝[37.12524.4816.32]^T

按照技术方案中的公式14，计算得到k_G的偏度系数λ_{X2_3}：

λ_{X2_3}＝[-2.6667-2.6667-2.6667]

按照技术方案中的公式15，计算得到k_G的峰度系数λ_{X2_4}：

λ_{X2_4}＝[8.11118.11118.1111]

III)负荷的单位调节功率；

按照技术方案中的公式16，将负荷单位调节功率标幺值区间[0,3]转换为各节点有名值区间(单位MW/Hz)：

$\left[\begin{array}{cc}{a}_5& b_5\\ a_7& b_7\\ a_9& b_9\end{array}\right]=\left[\begin{array}{cc}0& 5.4\\ 0& 6\\ 0& 7.5\end{array}\right]$

按照技术方案中的公式17，计算得到负荷的单位调节功率k_L的期望μ_X3(单位MW/Hz)，

μ_X3＝[2.733.75]^T

按照技术方案中的公式18，计算得到k_L的均方差σ_X3：

σ_X3＝[1.55881.7322.165]^T

按照技术方案中的公式19，计算得到k_L的偏度系数λ_{X3_3}：

λ_{X3_3}＝[000]^T

按照技术方案中的公式20，计算得到k_L的峰度系数λ_{X3_4}：

λ_{X3_4}＝[9/59/59/5]^T

5)根据2m+1点估计法构造随机变量的样本点，确定样本值及其权重系数；

完成步骤4后，根据步骤4得到的统计特征值，计算每个随机变量的样本值和对应的权重系数。

先按照技术方案中的公式22计算随机变量位置系数和权重系数，然后根据公式23计算随机变量的三个样本值。

a)负荷扰动；

计算得到扰动预测误差ΔP_{L_pre}的三个位置系数：

${\mathrm{\ξ}}_{X1}=\left[\begin{array}{ccc}1.7132& -1.7132& 0\\ 1.7132& -1.7132& 0\\ 1.7132& -1.7132& 0\end{array}\right]$

计算得到预测误差ΔP_{L_pre}的三个权重系数

$w_{X1}=\left[\begin{array}{ccc}0.1667& 0.1667& -0.2222\\ 0.1667& 0.1667& -0.2222\\ 0.1667& 0.1667& -0.2222\end{array}\right]$

按照技术方案中的公式21计算得到预测误差ΔP_{L_pre}的三个样本值

${Z^{\′}}_{X1}=\left[\begin{array}{ccc}1.5588& -1.5588& 0\\ 1.7321& -1.7321& 0\\ 2.1651& -2.1651& 0\end{array}\right]$

其中Z'_X1(i,j)表示第i个功率扰动预测误差的第j个样本值，i＝1,2,3，对应节点序号5,7,9；j＝1,2,3。

接着，按照技术方案中的公式7计算得到负荷扰动的样本值

$Z_{X1}=\left[\begin{array}{ccc}10.5588& 7.4412& 9\\ 11.7321& 8.2679& 10\\ 14.6651& 10.3349& 12.5\end{array}\right]$

其中Z_X1(i,j)表示第i个功率扰动的第j个样本值，i＝1,2,3，对应节点序号5,7,9；j＝1,2,3。

b)发电机的单位调节功率；

计算得到发电机的单位调节功率k_G的三个位置系数：

${\mathrm{\ξ}}_{X2}=\left[\begin{array}{ccc}0.3333& -3& 0\\ 0.3333& -3& 0\\ 0.3333& -3& 0\end{array}\right]$

计算得到k_G的三个权重系数

$w_{X2}=\left[\begin{array}{ccc}0.9000& 0.1& -0.8889\\ 0.9000& 0.1& -0.8889\\ 0.9000& 0.1& -0.8889\end{array}\right]$

计算得到k_G的三个样本值

$Z_{X2}=\left[\begin{array}{ccc}123.75& 0& 111.375\\ 81.6& 0& 73.44\\ 54.4& 0& 48.96\end{array}\right]$

其中Z_X2(i,j)表示第i个发电机的单位调节功率的第j个样本值，i＝1,2,3，j＝1,2,3。

c)负荷的单位调节功率；

计算得到负荷的单位调节功率k_L的三个位置系数：

${\mathrm{\ξ}}_{X3}=\left[\begin{array}{ccc}1.3416& -1.3416& 0\\ 1.3416& -1.3416& 0\\ 1.3416& -1.3416& 0\end{array}\right]$

计算得到k_L的三个权重系数

$w_{X3}=\left[\begin{array}{ccc}0.2778& 0.2778& -0.4444\\ 0.2778& 0.2778& -0.4444\\ 0.2778& 0.2778& -0.4444\end{array}\right]$

计算得到k_L的三个样本值

$Z_{X3}=\left[\begin{array}{ccc}4.7914& 0.6086& 2.7\\ 5.3238& 0.6762& 3\\ 6.6547& 0.8453& 3.75\end{array}\right]$

其中Z_X3(i,j)表示第i个负荷的单位调节功率的第j个样本值，

i＝1,2,3，对应节点序号5,7,9；j＝1,2,3。

上述计算完成后，根据技术方案中的步骤5中的系统样本构造方法，进一步得到系统的3m个样本.

6)参数初始化；

完成步骤5后，设置初始参数。最大迭代次数Tmax＝20，收敛精度ε＝10^-5。

按照技术方案中的公式24，计算得到基准状态下节点的注入无功功率(单位：MVar)

$Q_{input}^0={\left[\begin{array}{ccccccccc}24.0690& 14.4601& 3.6490& 0& -30& 0& -35& 0& -50\end{array}\right]}^T$

7)采用牛拉法计算各样本对应的系统状态变量值；

完成步骤6后，分别计算3m个样本点对应的确定性潮流，获得该样本点对应的系统频率f，节点电压幅值U和相角θ，发电机的实际有功出力P_G和负荷吸收的实际有功功率P_L。

下面以样本Z(1,3)为例，给出步骤7的相关中间结果。其中前三个数分别为节点5,7,9上的负荷扰动值，中间三个数为发电机的单位调节功率，最后三个数分别为节点5,7,9上负荷的单位调节功率。

i)计算初始不平衡量；

首先，设置迭代次数初值，令t＝0。

接着，根据样本，得到各节点的实际负荷扰动(单位：MW)：

ΔP_L＝[00009010012.5]^T

然后，按照技术方案中的公式27，计算此状态下个节点的注入有功功率(单位：MW)：

$P_{i\_input}^0=\[\begin{array}{ccccccccc}71.9547& 163& 85& 0& -99& 0& -110& 0& -137.5\end{array}\]$

接着，按照技术方案中公式25计算得到初始不平衡量。

$\left\{\begin{array}{c}{\mathrm{\ΔP}}^0=\[\begin{array}{ccccccccc}0& 0& 0& 0& -0.09& 0& -0.1& 0& -0.125\end{array}\]\\ {\mathrm{\ΔQ}}^0=\[\begin{array}{ccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\]\end{array}\right.$

ii)计算雅可比矩阵；

结合步骤1输入的数据，计算初次迭代的雅可比矩阵，按照技术方案中的公式29-38。

迭代次数自增，t＝t+1＝1。

iii)解修正方程，计算修正量；

步骤ii完成后，计算待求修正后的新值，公式按照技术方案中的公式39-42。

步骤ii完成后，按照技术方案中的公式39，计算得到第t＝1次迭代的修正量ΔX

其中第一部分是电压相角修正量(单位：弧度)，共8个；第二部分是与电压幅度相关修正量(标幺值)，共6个；第三部分为系统频率(标幺值)，共一个。

完成上述计算后，按照技术方案中的公式41和42，计算变量新值：

计算得到电压幅值新值(标幺值)：

Um¹＝[1110.98570.97281.00240.98370.99480.9546]

计算得到电压相角新值(单位：弧度)：

θ¹＝[00.16560.0772-0.0506-0.08470.0237-0.00300.0569-0.0922]计算得到频率新值(单位Hz)：

f¹＝49.8680

iv)计算不平衡量并进行收敛判断；

完成步骤iii后，计算新值状态下的不平衡量并判断是否收敛，公式按照技术方案中的公式25-28。

完成步骤iii后，按照技术方案中的公式25-28，计算得到修正后新值对应的节点功率的不平衡量ΔF：

完成上述计算后，按照技术方案中的公式44判断是否收敛：

max{|ΔF|}＝0.0012＞ε

收敛不满足，转入到步骤ii进行下一次迭代.

v)结果输出；

经过t＝3次迭代后，收敛完成。得到该系统样本点对应的潮流分布信息。

系统频率f(单位Hz)：

f＝49.8679

节点电压幅值Um(标幺值)：

Um＝[1110.98560.97271.00240.98370.99470.9545]

节点电压相角θ(单位：度)：

θ＝[09.48984.4198-2.9035-4.85901.3546-0.17433.2604-5.2883]

节点的实际有功负荷P_L(单位MW)：

P_L＝[000098.64320109.60360137.0045]

发电机的实际有功出力P_G(单位：MW)：

P_G＝[86.6710172.703891.4692]

8)根据所得结果及其对应权重系数，计算变量的统计特征值；

完成系统所有样本潮流计算之后，按照技术方案中的公式45-47，计算得到变量的期望。

系统频率期望f：

f＝49.8596(Hz)

系统状态变量如表I：

表I.扰动后的系统状态

实验效果

为验证本方法的有效性和正确性，采用蒙特卡洛方法计算上述算例，并将结果与本方法结果相对比。同时，为了验证计及系统一次调频作用及其系数不确定性的必要性，对本算例分别计算不计一次调频作用、计及一次调频作用但系数确定两种方案的系统状态数据，并与本发明所提方法的数据进行对比。

1)与蒙特卡洛方法对比

采用蒙特卡洛方法，抽样10000次，计算本算例，得到系统状态数据。并与本发明计算方法所得数据进行比较，按照下式计算相对误差ε，验证其有效性和正确性。

$\mathrm{\ϵ}=\frac{|x_{MCS}-x_{PEM}|}{\left|x_{MCS}\right|}\×100\%$

其中x_MCS和x_PEM分别代表蒙特卡洛计算数据和本方法计算数据。表II给出了两种方法的计算出的节点电压幅值和相角，以及两者的相对误差。

表II.两种方法计算的节点电压及误差分析

由表II可以看出，本方法所得结果与蒙特卡洛抽样10000次的结果最大相对误差为1.5％，最小为0。由此证明了本方法的有效性和正确性。

2)与常规潮流模型和系数确定的动态潮流模型计算结果对比

为验证本发明的效果，给出了方案A和方案B的计算结果，并与本发明所提方案进行对比。其中方案C是本专利所提出的计及一次调频不确定性的概率潮流分析方法。方案A是常规概率潮流模型，即不考虑系统的一次调频特性，认为系统频率固定为50Hz。方案B是动态潮流模型，考虑系统一次调频特性，且认为系统调频参数为确定值。发电机的单位调节功率标幺值取25，负荷的单位调节功率标幺值取1.5，两种方案中的负荷扰动以及其余参数与上述算例完全一样。

表III给出了方案A、方案B和方案C计算出的系统频率，表IV给出了方案A、方案B和方案C计算出的发电机有功出力，表V给出了方案A、方案B和方案C计算出的各节点的有功负荷。

表III.系统频率f(单位：Hz)

系统频率	方案A	方案B	方案C
				f	50	49.8806	49.8596

表IV.发电机有功出力PG(单位：MW)

节点序号

方案A

方案B

方案C

1	103.4311	86.7296	86.2190
				2	163	172.7425	172.8777
3	85.	91.4950	91.7016

表V.节点有功负荷PL(单位：MW)

节点序号	方案A	方案B	方案C
				5	99	98.6776	98.6222
7	110	109.6418	109.5803
				9	137.5	137.0523	136.9760

由上表可以看出，是否计及系统的一次调频特性对潮流计算结果有很大的影响。特别对于系统原平衡节点(节点1)，计及一次调频特性得到的发电机有功出力与常规概率潮流计算得到的发电机有功出力存在很大的差别。同时还可看出，发电机和负荷单位调节功率的不确定性对潮流计算结果影响明显。计及系统一次调频特性及其系数不确定性的概率潮流方法所得潮流结果可以更精确地描述扰动时系统的潮流分布。

Claims (1)

1.一种计及一次调频不确定性的概率潮流分析方法，其特征在于，包括以下步骤：

1)输入确定参数及随机变量的分布信息；

输入常规潮流计算所需要的基础信息、随机变量的分布信息和随机变量的数量；

所述输入常规潮流计算所需要的基础信息，包括：网络拓扑数据，支路阻抗，母线并联电导电纳，节点类型，节点负荷功率，PV节点上发电机出力数据，PQ节点上发电机出力数据，发电机的额定容量及功率因素，发电机的调差系数标幺值R_g和系统负荷扰动功率预测值；

所述随机变量的分布信息，随机变量包括：节点负荷扰动预测误差ΔP_{L_pre}、负荷的单位调节功率k_L和发电机的单位调节功率k_G；所述节点负荷扰动预测误差ΔP_{L_pre}服从正态分布，确定其正态分布的期望和方差；所述发电机的单位调节功率k_G服从二项分布，确定发电机单位调节功率标幺值的两点取值K_Gg＝1/R_g和K_Gg＝0，及其对应的概率；所述负荷的单位调节功率k_L服从均匀分布，确定其标幺值k_Lg的分布区间[a',b']；

所述随机变量的数量为输入系统中随机变量的数量，总个数m，其中随机变量ΔP_{L_pre}的个数为m_p,随机变量k_L的个数为m_kL，随机变量k_G的个数为m_kG；

2)计算基准状态下的潮流分布；

结合步骤1获得的所述常规潮流基础信息，建立如下常规潮流模型；

$\left\{\begin{array}{c}{P}_{Gi}-P_{Li}^0=U_i\underset{j=1}{\overset{j=n}{\Σ}}{U}_j(G_{ij}{\mathrm{cos\δ}}_{ij}+B_{ij}{\mathrm{sin\δ}}_{ij})\\ Q_{Gi}-Q_{Li}^0=U_i\underset{j=1}{\overset{j=n}{\Σ}}{U}_j(G_{ij}{\mathrm{sin\δ}}_{ij}-B_{ij}{\mathrm{cos\δ}}_{ij})\end{array}\right.---\left(1\right)$

公式1是基于极坐标的常规潮流方程，其中P_Gi,Q_Gi分别表示节点i所接发电机的有功出力和无功出力；分别表示基准状态下，节点i上负荷吸收的有功功率和无功功率；U_i,U_j分别表示节点i的电压幅值和节点j的电压幅值，δ_ij＝θ_i-θ_j表示节点i、j间电压的相角差，其中θ_i和θ_j分别表示节点i和节点j电压的相角；G_ij,B_ij分别表示节点i、j间的电导和电纳，若节点i、j间无支路则对应的导纳为零；

运用牛拉法求解上述方程，得到基准状态下发电机有功出力无功出力节点电压幅值和相角

3)构建计及系统一次调频特性的动态潮流模型；

结合步骤1获得的所述常规潮流基础信息以及步骤2求解得到的所述基准状态变量，构建负荷扰动时，计及负荷和发电机一次调频作用的动态潮流模型如下：

$\left\{\begin{array}{c}{P}_{Gi}-P_{Li}=U_i\underset{j=1}{\overset{j=n}{\Σ}}{U}_j(G_{ij}{\mathrm{cos\δ}}_{ij}+B_{ij}{\mathrm{sin\δ}}_{ij})\\ Q_{Gi}^0-Q_{Li}^0=U_i\underset{j=1}{\overset{j=n}{\Σ}}{U}_j(G_{ij}{\mathrm{sin\δ}}_{ij}-B_{ij}{\mathrm{cos\δ}}_{ij})\\ P_{Gi}=P_{Gi}^0-k_{Gi}\mathrm{\Δ}f\\ P_{Li}=P_{Li}^{\′}+k_{Li}\mathrm{\Δ}f=(P_{Li}^0+{\mathrm{\ΔP}}_{Li})+k_{Li}\mathrm{\Δ}f\end{array}\right.---\left(2\right)$

公式2中，系统发生扰动后，计及一次调频作用时节点i上发电机的实际有功出力用P_Gi表示，负荷吸收的实际有功功率用P_Li表示；k_Gi表示节点i上发电机的单位调节功率,k_Li表示负荷的单位调节功率；Δf＝f-f₀为系统的频率偏移，其中f₀和f分别表示基准状态的系统频率(50Hz)和扰动后的系统频率；P′_Li表示扰动发生瞬间节点i上的有功，表示扰动前基准状态节点i上的有功负荷,ΔP_Li表示基准状态节点i上发生的实际负荷扰动功率；

4)根据分布信息，计算随机变量所需的统计特征值；

对概率密度为f_X(x)的变量X,求取下列统计特征值：

期望μ_X：

${\mathrm{\μ}}_X=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{\mathrm{xf}}_X\left(x\right)dx---\left(3\right)$

标准差σ_X：

${\mathrm{\σ}}_X=\sqrt{\left(\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{(x-{\mathrm{\μ}}_X)}^2{f}_X\left(x\right)dx\right)}---\left(4\right)$

X的j阶中心距M_j(X)：

$M_j\left(X\right)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}{(x-{\mathrm{\μ}}_X)}^j{f}_X\left(x\right)dx---\left(5\right)$

定义λ_{X_j}表示随机变量X的j阶中心距与标准差σ_X的j次方的比值：

${\mathrm{\λ}}_{X\_j}=\frac{M_j\left(X\right)}{{\mathrm{\σ}}_X^j}---\left(6\right)$

其中，λ_{X_1}＝0，λ_{X_2}＝1，λ_{X_3}和λ_{X_4}分别表示随机变量X的偏度系数和峰度系数；

结合公式3-6，并根据步骤1获得的所述随机变量分布信息，计算所述随机变量的期望、均方差、偏度系数和峰度系数；计算过程如下：

I)负荷扰动预测误差；

实际的扰动功率ΔP_L由负荷扰动的预测值P_{L_pre}及预测误差ΔP_{L_pre}两部分组成，计算表达式如公式7：

ΔP_L＝P_{L_pre}+ΔP_{L_pre}(7)

用变量X1表示预测误差ΔP_{L_pre}，各统计特征值求取如下：

利用公式8计算扰动预测误差ΔP_{L_pre}的偏度系数λ_{X1_3}：

${\mathrm{\λ}}_{X1\_3}=\frac{M_3\left(X1\right)}{{\mathrm{\σ}}_{X1}^3}=\frac{0}{{\mathrm{\σ}}_{X1}^3}=0---\left(8\right)$

利用公式9计算扰动预测误差ΔP_{L_pre}的峰度系数λ_{X1_4}：

${\mathrm{\λ}}_{X1\_4}=\frac{M_4\left(X1\right)}{{\mathrm{\σ}}_{X1}^4}=\frac{3{\mathrm{\σ}}_{X1}^4}{{\mathrm{\σ}}_{X1}^4}=3---\left(9\right)$

II)发电机的单位调节功率；

机组的单位调节功率k_G服从二项分布，p表示参与一次调频的概率，对应的单位调节功率的值为1/R，其中R是机组的调差系数；q＝1-p表示不参与一次调频的概率，对应的单位调节功率的值为0；

将步骤1中输入的调差系数标幺值R_g转换为其有名值R，转换公式如下：

$R=R_g\frac{f_N}{P_{GN}}---\left(10\right)$

其中f_N＝50Hz，为系统额定频率；P_GN是机组的额定有功功率，通过公式11求取：

P_GN＝S_GNcosα(11)

其中S_GN机组的额定容量，α为功率因素；

用变量X2表示发电机的单位调节功率k_G，其各统计特征值求取如下：

利用公式12计算发电机的单位调节功率k_G的期望μ_X2：

${\mathrm{\μ}}_{X2}=p\frac{1}{R}+q\×0=\frac{p}{R}---\left(12\right)$

利用公式13计算发电机的单位调节效率k_G的均方差σ_X2：

${\mathrm{\σ}}_{X2}=\sqrt{p{(\frac{1}{R}-\frac{p}{R})}^2+q{(0-\frac{p}{R})}^2}=\frac{\sqrt{pq}}{R}---\left(13\right)$

利用公式14计算发电机的单位调节功率k_G的偏度系数λ_{X2_3}：

${\mathrm{\λ}}_{X2\_3}=\frac{M_3\left(X2\right)}{{\mathrm{\σ}}_{X2}^3}=\frac{p{(1/R-p/R)}^3+q{(0-p/R)}^3}{{(\sqrt{\mathrm{pq}}/R)}^3}=\frac{q^2-p^2}{\sqrt{pq}}---\left(14\right)$

利用公式15计算发电机的单位调节功率k_G的峰度系数λ_{X2_4}：

${\mathrm{\λ}}_{X2\_4}=\frac{M_4\left(X2\right)}{{\mathrm{\σ}}_{X2}^4}=\frac{p{(1/R-p/R)}^4+q{(0-p/R)}^4}{{(\sqrt{pq}/R)}^4}=\frac{q^3+p^3}{pq}---\left(15\right)$

III)负荷的单位调节功率；

负荷的单位调节功率标幺值k_Lg在[a',b']区间内均匀分布；与步骤II一样，需要先将标幺值区间[a',b']转换为有名值区间[a,b]，按照公式16转换：

$k_L=k_{Lg}\frac{P_{LN}}{f_N}---\left(16\right)$

其中P_LN表示额定有功；f_N＝50Hz，为系统额定频率；

用变量X3表示负荷的单位调节功率k_L，其各统计特征值求取如下：

利用公式17计算负荷的单位调节功率k_L的期望μ_X3：

${\mathrm{\μ}}_{X3}=\frac{b-a}{2}---\left(17\right)$

利用公式18计算负荷的单位调节功率k_L的均方差σ_X3：

${\mathrm{\σ}}_{X3}=\sqrt{\underset{a}{\overset{b}{\∫}}{(x-{\mathrm{\μ}}_{X3})}^2\frac{1}{b-a}dx}=\frac{b-a}{2\sqrt{3}}---\left(18\right)$

利用公式19计算负荷的单位调节功率k_L的偏度系数λ_{X3_3}：

${\mathrm{\λ}}_{X3\_3}=\frac{M_3\left(X3\right)}{{\mathrm{\σ}}_{X3}^3}=\frac{\underset{a}{\overset{b}{\∫}}{(x-{\mathrm{\μ}}_{X3})}^3\frac{1}{b-a}dx}{{\mathrm{\σ}}_{X3}^3}=0---\left(19\right)$

利用公式20计算负荷的单位调节功率k_L的峰度系数λ_{X3_4}：

${\mathrm{\λ}}_{X3\_4}=\frac{M_4\left(X3\right)}{{\mathrm{\σ}}_{X3}^4}=\frac{\underset{a}{\overset{b}{\∫}}{(x-{\mathrm{\μ}}_{X3})}^4\frac{1}{b-a}dx}{{\mathrm{\σ}}_{X3}^4}=\frac{9}{5}---\left(20\right)$

5)根据2m+1点估计法构造随机变量的样本点，确定样本值及其权重系数；

首先求取各样本值对应的权重系数；根据点估计原理，对系统中的所有随机变量z_i(z_i∈{z₁,z₂,…z_m})分别按照公式21构造三个样本值z_i,1,z_i,2,z_i,3；

z_i,j＝μ_i+ξ_i,jσ_ii＝1,…m；j＝1,2,3(21)

其中ξ_i,j表示第i个随机变量的第j个位置系数，通过公式22求取；

$\left\{\begin{array}{c}{\mathrm{\ξ}}_{i,1}=\frac{{\mathrm{\λ}}_{i,3}}{2}+\sqrt{{\mathrm{\λ}}_{i,4}-\frac{3}{4}{\mathrm{\λ}}_{i,3}^2}\\ {\mathrm{\ξ}}_{i,2}=\frac{{\mathrm{\λ}}_{i,3}}{2}-\sqrt{{\mathrm{\λ}}_{i,4}-\frac{3}{4}{\mathrm{\λ}}_{i,3}^2}\\ {\mathrm{\ξ}}_{i,3}=0\end{array}\right.---\left(22\right)$

其中λ_i3和λ_i4分别表示步骤4中求取的第i个随机变量的偏度系数和峰度系数；

基于公式22得到的参数，利用公式23求取各样本值对应的权重系数w_i,1、w_i,2和w_i,3；

$\left\{\begin{array}{c}{w}_{i,1}=\frac{1}{{\mathrm{\ξ}}_{i1}^2-{\mathrm{\ξ}}_{i1}{\mathrm{\ξ}}_{i2}}\\ w_{i,2}=\frac{1}{{\mathrm{\ξ}}_{i2}^2-{\mathrm{\ξ}}_{i1}{\mathrm{\ξ}}_{i2}}\\ w_{i,3}=\frac{1}{m}-\frac{1}{{\mathrm{\λ}}_{i,4}-{\mathrm{\λ}}_{i,3}^2}\end{array}\right.---\left(23\right)$

其中m表示系统中随机变量的个数；

然后，基于得到的负荷扰动预测误差的三个样本点，按照步骤4中的公式7，计算得到负荷扰动的样本；

最后构造系统样本点；构造随机变量样本值z_i,j(i＝1,…m；j＝1,2,3)对应的系统样本点时，随机变量z_i＝z_i,j，其余变量取均值，即样本点Z(i,j)＝(μ₁,…,z_i,j,…μ_m)；其中Z(1,3)＝Z(2,3)＝…＝Z(m,3)＝(μ₁,…,μ_i,…μ_m)，只需进行一次计算，故对有m个随机变量的系统，只需完成2m+1次确定性潮流计算即可；

6)参数初始化；

完成步骤5后，确定潮流计算所需的初始参数；设置最大迭代次数Tmax，收敛精度ε＝10^-5～10^-3；根据步骤2得到的基准状态数据，确定发电机有功出力初值、无功出力初值分别为和节点i电压幅值和相角初值分别为和系统频率初值f_N＝50Hz；

利用公式24计算节点i的注入无功

$Q_{i\_input}^0=Q_{Gi}^0-Q_{Li}^0---\left(24\right)$

7)采用牛拉法计算各样本对应的系统状态变量值；

完成步骤6后，基于步骤5得到的由负荷扰动、发电机单位调节功率和负荷单位调节功率组成的系统样本和步骤3构造的动态潮流模型，分别计算3m个样本点对应的确定性潮流，获得该样本点对应的系统频率f，节点电压幅值U和相角θ，发电机的实际有功出力P_G和负荷吸收的实际有功功率P_L，计算过程包括；

i)设置迭代次数初值，计算初始不平衡量；

首先设置迭代次数初值t＝0；

利用公式25计算基准状态发生扰动ΔP_Li瞬间的初始不平衡量；

$\left\{\begin{array}{cc}{\mathrm{\ΔP}}_i^0=P_{i\_input}^0-U_i^0\underset{j=1}{\overset{j=n}{\mathrm{\Σ}}}{U}_j^0(G_{ij}{\mathrm{cos\δ}}_{ij}^0+B_{ij}{\mathrm{sin\δ}}_{ij}^0)& (i\∈\[pq;pv;ref\])\\ {\mathrm{\ΔQ}}_i^0=Q_{i\_input}^0-U_i^0\underset{j=1}{\overset{j=n}{\mathrm{\Σ}}}{U}_j^0(G_{ij}{\mathrm{sin\δ}}_{ij}^0-B_{ij}{\mathrm{cos\δ}}_{ij}^0)& (i\∈pq)\end{array}\right.---\left(25\right)$

公式25中为基准状态节点i和节点j电压的相角差，即；

${\mathrm{\δ}}_{ij}^0={\mathrm{\θ}}_i^0-{\mathrm{\θ}}_j^0---\left(26\right)$

表示扰动后节点i的注入有功，根据公式27求取；

$\{\begin{array}{c}{P}_{i\_input}^0=P_{Gi}^0-{P^{\′}}_{Li}^0\\ {P^{\′}}_{Li}^0=P_{Li}^0+{\mathrm{\ΔP}}_{Li}=P_{Li}^0+P_{Li\_pre}+{\mathrm{\ΔP}}_{Li\_pre}\end{array}---\left(27\right)$

其中和分别表示基准状态下节点i上发电机的有功出力和有功负荷；表示扰动瞬间节点i的有功负荷；P_{Li_pre}和ΔP_{Li_pre}分别表示扰动预测值和扰动预测误差；

此时系统的不平衡量为ΔF；

ii)计算雅克比矩阵；

基于步骤3提出的动态潮流模型，计算雅克比矩阵；

$J=\left[\begin{array}{ccc}H& N& M\\ K& L& T\\ W& O& V\end{array}\right]---\left(29\right)$

其中H为n_pq+n_pv阶的方阵；

$H_{ij}=\frac{\∂{\mathrm{\ΔP}}_i}{\∂{\mathrm{\θ}}_j}=\left\{\begin{array}{cc}{U}_i^t{U}_j^t(G_{ij}{\mathrm{sin\δ}}_{ij}^t-B_{ij}{\mathrm{cos\δ}}_{ij}^t)& (i\≠j)\\ -U_i^t\underset{k\≠i}{\overset{n}{\underset{k=1}{\mathrm{\Σ}}}}{U}_k^t(G_{ik}{\mathrm{sin\δ}}_{ik}^t-B_{ik}{\mathrm{cos\δ}}_{ik}^t)& (i=j)\end{array}\right.,(i\∈\[pq;pv\];j\∈\[pq;pv\])---\left(30\right)$

N为(n_pq+n_pv)×n_pq阶的矩阵；

$H_{ij}=\frac{\∂{\mathrm{\ΔP}}_i}{\∂U_j}{U}_j=\left\{\begin{array}{cc}{U}_i^t{U}_j^t(G_{ij}{\mathrm{cos\δ}}_{ij}^t+B_{ij}{\mathrm{sin\δ}}_{ij}^t)& (i\≠j)\\ U_i^t\underset{k=1}{\overset{n}{\Σ}}{U}_k^t(G_{ik}{\mathrm{cos\δ}}_{ik}^t+B_{ik}{\mathrm{sin\δ}}_{ik}^t)+{\left(U_i^t\right)}^2{G}_{ii}& (i=j)\end{array}\right.,(i\∈\[pq;pv\];j\∈pq)---\left(31\right)$

M为(n_pq+n_pv)×1阶的矩阵；

$M_i=\frac{\∂{\mathrm{\ΔP}}_i}{\∂f}=(k_{Gi}+k_{Li}),(i\∈\[pq:pv\])---\left(32\right)$

K为n_pq×(n_pq+n_pv)阶的矩阵；

$H_{ij}=\frac{\∂{\mathrm{\ΔQ}}_i}{\∂{\mathrm{\θ}}_j}=\left\{\begin{array}{cc}{-U}_i^t{U}_j^t(G_{ij}{\cos \mathrm{\δ}}_{ij}^t+B_{ij}{\sin \mathrm{\δ}}_{ij}^t)& (i\≠j)\\ U_i^t\underset{k\≠i}{\overset{n}{\underset{k=1}{\mathrm{\Σ}}}}{U}_k^t(G_{ik}{\cos \mathrm{\δ}}_{ik}^t+B_{ik}{\sin \mathrm{\δ}}_{ik}^t)& (i=j)\end{array}\right.,(i\∈pq;j\∈\[pq;pv\])---\left(33\right)$

L为n_pq×n_pq阶的矩阵；

$L_{ij}=\frac{\∂{\mathrm{\ΔQ}}_i}{\∂U_j}{U}_j=\left\{\begin{array}{cc}{U}_i^t{U}_j^t(G_{ij}{\mathrm{sin\δ}}_{ij}^t-B_{ij}{\mathrm{cos\δ}}_{ij}^t)& (i\≠j)\\ U_i^t\underset{k=1}{\overset{n}{\Σ}}{U}_k^t(G_{ik}{\mathrm{sin\δ}}_{ik}^t-B_{ik}{\mathrm{cos\δ}}_{ik}^t)-{\left(U_i^t\right)}^2{B}_{ii}& (i=j)\end{array}\right.,(i\∈pq;j\∈pq)---\left(34\right)$

T为n_pq×1阶的矩阵；

$T_i=\frac{\∂Q_i}{\∂f}=0,(i\∈pq)---\left(35\right)$

W为1×(n_pq+n_pv)阶的矩阵；

$W_j=\frac{\∂{\mathrm{\ΔP}}_i}{\∂{\mathrm{\θ}}_j}=U_i^t{U}_j^t(G_{ij}{\mathrm{sin\δ}}_{ij}^t-B_{ij}{\mathrm{cos\δ}}_{ij}^t),(i=ref,j\∈\[pq;pv\])---\left(36\right)$

O为1×n_pq阶的矩阵；

$O_j=\frac{\∂{\mathrm{\ΔP}}_i}{\∂U_j}{U}_j=U_i^t{U}_j^t(G_{ij}{\mathrm{cos\δ}}_{ij}^t+B_{ij}{\mathrm{sin\δ}}_{ij}^t),(i=ref;j\∈pq)---\left(37\right)$

v为1×1阶的矩阵；

$V=\frac{\∂{\mathrm{\ΔP}}_{ref}}{f}=k_{Gi}+k_{Li},(i=ref)---\left(38\right)$

其pq,pv分别表示PQ节点序号矩阵，PV节点序号矩阵；ref为平衡节点的序号；n,n_pv,n_pq,n_ref(n_ref＝1)分别表示系统总的节点个数，PQ节点个数，PV节点个数和平衡节点个数；

迭代次数自增，即t＝t+1；

iii)解修正方程，计算修正量；

基于步骤i和步骤ii，求解修正方程39，得到修正量；

ΔX＝-J^-1ΔF(39)

ΔX为变量的修正量，参见公式40；

利用公式41对电压修正量进行变换；

${\mathrm{\ΔU}}_i=\frac{{\mathrm{\ΔU}}_i}{U_i}{U}_i^{t-1},(i\∈pq)---\left(41\right)$

其中表示上次迭代得到的节点i的电压幅值，初次迭代时指给定的电压幅值初值；

进一步，结合公式40和41，计算得到的修正量修正原有的变量，得到新值；

X^t＝X^t-1+ΔX(42)

iv)计算不平衡量并进行收敛判断；

按照公式25计算新值下的不平衡量ΔF；

判断不平衡量ΔF是否满足收敛判据；

max{|ΔF|}≤ε(44)

若收敛判据满足则迭代结束转到步骤v；若收敛判据不满足且迭代次数t＜Tmax则跳转至步骤ii，继续迭代直至满足收敛条件；若判据不满足且t＝Tmax，说明迭代不收敛，跳出循环，转至步骤v；

v)结果输出；

判断此时t是否满足t＜Tmax，若满足则此时得到的结果X即为样本Z(i,j)对应的状态变量值；否则表示潮流迭代不收敛；

8)计算变量x的统计特征值；

基于步骤7所得的3m个系统样本计算出的潮流结果以及步骤5计算得到的各样本点对应的权重系数，根据公式45-47计算变量的期望和方差；

$E\[x^k\]=\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{3}{\Σ}}{w}_{Lj}{\left(X\right(i,j\left)\right)}^k---\left(45\right)$

其中

X(1,3)＝X(2,3)＝…＝X(m,3)(46)

$\left\{\begin{array}{c}{\mathrm{\μ}}_x=E\[x\]\\ {\mathrm{\σ}}_x=\sqrt{E\[x^2\]-{{\mathrm{\μ}}_x}^2}\end{array}\right.---\left(47\right)$

其中X(i,j)表示样本Z(i,j)计算得到的变量X的值；μ_x和σ_x分别表示变量X的期望和均方差；变量X包含节点电压相角，电压幅值和频率，