[发明专利]一种计及一次调频不确定性的概率潮流分析方法

权利要求书

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获得的所述常规潮流基础信息，建立如下常规潮流模型；

PGi-PLi0=UiΣj=1j=nUj(Gijcosδij+Bijsinδij)QGi-QLi0=UiΣj=1j=nUj(Gijsinδij-Bijcosδij)---(1)]]>

公式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求解得到的所述基准状态变量，构建负荷扰动时，计及负荷和发电机一次调频作用的动态潮流模型如下：

PGi-PLi=UiΣj=1j=nUj(Gijcosδij+Bijsinδij)QGi0-QLi0=UiΣj=1j=nUj(Gijsinδij-Bijcosδij)PGi=PGi0-kGiΔfPLi=PLi′+kLiΔf=(PLi0+ΔPLi)+kLiΔf---(2)]]>

公式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：

μX=∫-∞+∞xfX(x)dx---(3)]]>

标准差σ_X：

σX=(∫-∞+∞(x-μX)2fX(x)dx)---(4)]]>

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

Mj(X)=∫-∞+∞(x-μX)jfX(x)dx---(5)]]>

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

λX_j=Mj(X)σXj---(6)]]>

其中，λ_{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}：

λX1_3=M3(X1)σX13=0σX13=0---(8)]]>

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

λX1_4=M4(X1)σX14=3σX14σX14=3---(9)]]>

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

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

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

R=RgfNPGN---(10)]]> 2

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

P_GN＝S_GNcosα (11)

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

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

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

μX2=p1R+q×0=pR---(12)]]>

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

σX2=p(1R-pR)2+q(0-pR)2=pqR---(13)]]>

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

λX2_3=M3(X2)σX23=p(1/R-p/R)3+q(0-p/R)3(pq/R)3=q2-p2pq---(14)]]>

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

λX2_4=M4(X2)σX24=p(1/R-p/R)4+q(0-p/R)4(pq/R)4=q3+p3pq---(15)]]>

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

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

kL=kL·PLNfN---(16)]]>

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

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

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

μX3=b-a2---(17)]]>

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

σX3=∫ab(x-μX3)21b-adx=b-a23---(18)]]>

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

λX3_3=M3(X3)σX33=∫ab(x-μX3)31b-adxσX33=0---(19)]]>

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

λX3_4=M4(X3)σX34=∫ab(x-μX3)41b-adxσX34=95---(20)]]> 3

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σ_i i＝1,…m；j＝1,2,3 (21)

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

ξi,1=λi,32+λi,4-34λi,32ξi,2=λi,32-λi,4-34λi,32ξi,3=0---(22)]]>

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

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

wi,1=1ξi12-ξi1ξi2wi,2=1ξi12-ξi1ξi2wi,3=1m-1λi,4-λi,32---(23)]]>

其中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的注入无功

Qi_input0=QGi0-QLi0---(24)]]>

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

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

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

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

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

ΔPi0=Pi_input0-Ui0Σj=1j=nUj0(Gijcosδij0+Bijsinδij0)(i∈[pq;pv;ref])ΔQi0=Ui_input0-Ui0Σj=1j=nUj0(Gijsinδij0-Bijcosδij0)(i∈pq)---(25)]]>

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

δij0=θi0-θj0---(26)]]>

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

Pi_input0=PGi0-P′Li0P′Li0=PLi0+ΔPLi=PLi0+PLi_pre+ΔPLi_pre---(27)]]>

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

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

ii)计算雅克比矩阵；

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

J=HNMKLTWOV---(29)]]>

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

Hij=∂ΔPi∂θj=UitUjt(Gijsinδijt-Bijcosδijt)(i≠j)-UitΣk=1k≠inUkt(Giksinδikt-Bikcosδikt)(i=j)(i∈[pq;pv];j∈[pq;pv])---(30)]]>

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

Nij=∂ΔPi∂UjUj=UitUjt(Gijcosδijt+Bijsinδijt)(i≠j)UitΣk=1nUkt(Gikcosδikt+Biksinδikt)+(Uit)2Gii(i=j)(i∈[pq;pv];j∈pq)---(31)]]>

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

Mi=∂ΔPi∂f=(kGi+kLi)(i∈[pq:pv])---(32)]]>

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

Kij=∂ΔQi∂θj=-UitUjt(Gijcosδijt+Bijsinδijt)(i≠j)UitΣk=1k≠inUkt(Gikcosδikt+Biksinδikt)(i=j)(i∈pq;j∈[pq;pv])---(33)]]> 5

L为n_pq×n_pq阶的矩阵；

Lij=∂ΔQi∂UjUj=UitUjt(Gijsinδijt-Bijcosδijt)(i≠j)UitΣk=1nUkt(Giksinδikt-Bikcosδikt)-(Uit)2Bii(i=j)(i∈pq;j∈pq)---(34)]]>

T为n_pq×1阶的矩阵；

Ti=∂Qi∂f=0(i∈pq)---(35)]]>

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

Wj=∂ΔPi∂θj=UitUjt(Gijsinδijt-Bijcosδijt)(i=ref,j∈[pq;pv])---(36)]]>

O为1×n_pq阶的矩阵；

Oj=∂ΔPi∂UjUj=UitUjt(Gijcosδijt+Bijsinδijt)(i=ref;j∈pq)---(37)]]>

v为1×1阶的矩阵；

V=∂ΔPreff=kGi+kLi(i=ref)---(38)]]>

其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对电压修正量进行变换；

ΔUi=ΔUiUiUit-1(i∈pq)---(41)]]>

其中U_i^t-1表示上次迭代得到的节点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[xk]=Σi=1mΣj=13wi,j(X(i,j))k---(45)]]>

其中

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

μx=E[x]σx=E[x2]-μx2---(47)]]>

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